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1 . INTRODUCTION 

The use of ceramics in heavy duty diesel engine applications Is 
especially promising [1]. It has been shown in the TACOM/ Cummins Adiabatic 
Engine Program that reductions In heat losses at appropriate points in the 
diesel engine system result in substantially Increased exhaust enthalpy [2], 
One of the most promising concepts for taking advantage of this increased 
exhaust enthalpy is the turbocharged turbocoapounded diesel engine cycle [33. 

This engine concept consists of several sub-systems: compressor, 
reciprocator, turbocharger turbine, compounded turbine, ducting and heat 
exchangers. Opportunities for use of ceramic materials in the different sub- 
system components hare widely different degrees of difficulty and benefits. 
Therefore, there is a need for a computer simulation of this complete engine 
concept, at the appropriate level of detail, to enable engineers to define the 
trade-offs associated with Introducing ceramic materials in various parts of 
the total engine system, and to carry out system performance and optimization 
studies. 

This report describes the thermodynamic and heat transfer models used in 
a conputer program which simulates the behavior of the total engine system. 

The focus of this total system simulation i3 to define the mass flows and 
energy transfers (heat transfers, heat release, work transfers) in each sub- 
system and the relationship between the sub-systems. Slrsce this system model 
must function as a single-unit, a deliberate effort is made to maintain a 
balance in the complexity of the various sub-system descriptions. 

This effort has been funded by DOE over a two year period. A Phase I 
computer simulation of the turbocompounded engine system was developed during 
the first year. The mathematical models used in the Phase I simulation are 
described in NASA Contractors Report CR-17A755 [H]. In order to make the 
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Phase I program available as quickly as possible, relative simple heat 
transfer and conduction models were used for the major system components, A 
validation of this Phase I simulation against data from a turbocompounded 
engine system was successfully carried out and has also been reported [5], In 
the second phase of this program, a substantial extension and upgrading of 
critical sub-system models was carried out. Major areas of focus were: 
description of turbomachinery component performance; heat conduction through 
the piston, cylinder-head, cylinder-liner, exhaust port, exhaust manifold and 
ducting to incorporate material properties and wall designs appropriate to use 
of ceramics; radiation from the combustion gases; flow through the exhaust 
system. The present report describes the models used to represent the mass 
flow rates, heat transfers to and through the walls, work transfers, energy 
release In the combustion process, and thermodynamics of the working fluids 
properties and processes in this Phase II version of the turbocompounded 
diesel simulation. 

Figure 1 illustrates the overall model structure. The air flow is 
followed through an air filter, ducting, turbocharger compressor, ducting, 
cooler and engine intake system to the diesel reciprocator. The reciprocator 
simulation is a mathematical model of the processes occurring in the direct- 
injection four-stroke diesel-engine. Engine friction sub-models are then 
used to obtain brake quantities from the confuted reciprocator indicated 
performance quantities. The exhaust gas flow is followed through the engine 
exhaust ports, manifold runners and plenum, to the turbocharger turbine; it is 
then followed through additional ducting to the compounded turbine; the flow 
then passes through the exhaust system and muffler to the atmosphere. 

The report is arranged as follows. First, the basic assumptions of the 
reciprocator and the other component models are summarized. Then, the 
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mathematical relationships for mass and energy conservation for an open system 
(such as the manifolds and engine cylinder) are developed. Then, the models 
which represent each process in the individual components of the total system 
are described. Next, the method of solution of the complete model that 
results from the integration of the various sub-systems is presented, along 
with a typical set of program inputs and outputs. 




2. BASIC ASSUMPTIONS OF SYSTEM MODELS 
2.1 Reclprocator Engine Model 

The diesel four-stroke cycle is treated as a sequence of continuous 
processes: intake, compression, combustion (including expansion), and 

exhaust. The durations of the individual processes are as follows. The 
intake process begins then the intake valve opens (IVO) and ends when the 
Intake valve closes (IVC). The compression process begins at IVC and ends at 
the time of ignition (IGN). The combustion process begins when ignition 
occurs and ends When the exhaust valve opens (EVO). The exhaust process 
begins at EVO and ends at IVO (and not when the exhaust valve closes). 

In the reciprocator simulation, the system of interest is the 
instantaneous contents of a cylinder, i.e. air, fuel and combustion products. 
In general, this system is open to the transfer of mass, enthalpy, and energy 
in the form of work and heat. Throughout the cycle, the cylinder is treated 
as a variable volume plenum, spatially uniform in pressure. Furthermore, the 
cylinder contents are represented as one continuous medium by defining an 
average equivalence ratio and temperature in the cylinder at all times. 

Gas properties are calculated assuming ideal gas behavior. At low 
temperatures (below 1000 K), the cylinder contents are treated as a 
homogeneous mixture of non-reacting ideal gases. At high temperatures (above 
1000 K), the properties of the cylinder contents are calculated with allowance 
for chemical dissociation by assuming that the burned gases are in 
equilibrium, using an approximate calculation method based on hydrocarbon-air 
combustion. 

Quasi -steady, adiabatic, one-dimensional flow equations are used to 
predict mass flows past the intake and exhaust valves. The intake manifold 
and the exhaust port are treated as plenums whose pressure and temperature 
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history is determined by solution of the manifold state equations. When 
reverse flow past the intake valve occurs, rapid mixing of the back flow gases 
within the intake manifold is assumed. 

The compression process is defined so as to include the ignition delay 
period, i.e. the time interval between the start of the injection process (the 
point at wnich the injector needle starts to lift) and the ignition point (the 
start of positive heat release due to combustion). The total length of the 
ignition delay is related to the mean cylinder gas temperature and pressure 
during the delay period by an empirical Arrhenius expression. 

Combustion is modelled as a uniformly distributed heat release process. 

The rate of heat release is assumed to be proportional to the rate of fuel 
burning which is modelled empirically. Since the diesel combustion process is 
comprised of a pre-mixed and a diffusion-controlled combustion mechanism, 
Watson's fuel burning rate correlation [6], consisting of the sum of two 
algebraic functions, one for each combustion mechanism, is used. The fraction 
of the total fuel injected that is burnt by either mechanism depends on the 
length of the ignition delay period and the engine load and speed. 

Heat transfer is included in all the ermine processes. Convective heat 
transfer is modelled using available engine correlations based on turbulent flow 
in pipes. The characteristic velocity and length scales required to evaluate 
these correlations are obtained from a mean and turbulent kinetic energy model. 
Radiat'. se heat transfer is added during combustion. The steady-state inside 
vail surface temperatures of the piston, cylinder head, and liner can be either 
specified or calculated from a specification of the component vail structure. 
Additionally, the time-dependent temperature distribution in the piston and 
cylinder head can be computed using a one -dimensional finite difference model. 

Finally, an empirical friction model is used to convert the indicated 
engine performance quantities to brake performance quantities. 
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2.2 Other Component Models 

The reciprocator engine model calculates the state variables In one master 
cylinder of a multi-cylinder engine, while the manifolds and the other component 
models have inherent multi-cylinder capability. The Interaction between the 
master cylinder model and the other components is accounted for in the 
manifolds. To simulate the effect of additional cylinders on the manifold 
conditions, and hence on the entire system behavior, the conditions in the 
other cylinders are assumed to vary as echoes of the master cylinder, shifted by 
the appropriate phase angles. The following additional assumptions are made in 
order to develop models for the various components of the complete turbocharged 
turbocomp ounded diesel engine system: 

Intake air and exhaust g<is can be modelled as ideal gases. There is 
perfect and instantaneous mixing of all mass flows that enter each manifold (or 
section of the manifold) with the gases in the manifold. This implies that 
there is no spatial variation in properties within a manifold (or manifold 
section) at any instant of time, and all flows leaving a manifold have the 
properties of the manifold (or manifold section) contents. The connecting pipes 
to and from the Intake manifold are included as parts of the Intake manifold. 

The exhaust manifold is modeled as a series of connected open systems: exhaust 
port, manifold runner and manifold plenum. 

There is instantaneous mixing of all mass flows that enter the intake 
manifold with the gases in the manifold; thus there is no variation in 
properties wltnin the intake manifold at any Instant of time. All flows 
leaving the Intake manifold have the properties of the manifold contents. 
Similarly, the gas properties are assumed to be uniform within each individual 
section of the exhaust manifold at any Instant of time and the flow out of any 
section has the properties of that section. 




There are no mass transfers except along the routes Indicated; i.e., there 
are no losses or leakages from any component of the system. The flow through 
the wastegate can be specified as a fraction of the exhaust gas stream. The two 
. st gas streams, through and bypassing the turbocharger turbine, reunite 
with ro loss of thermal energy, at the turbocharger turbine exit pressure. An 
open-system duct model connects the turbocharger turbine and compounded turbine. 

The steady-state inside wall surface temperatures of the manifold sections 
and ducting can be specified or calculated from wall design information. 
Empirical correlations are used to calculate instantaneous heat transfer rates 
from the gas to the walls and pressure losses. Conduction through the walls is 
modeled as axisymraetric. The boundary conditions on the outside wall surface of 
each component are specified. From these, the steady-state temperature 
distribution within the walls of each component can be calculated. 

The turbomachinery performance is defined by maps that interrelate 
efficiency, pressure ratio, mass flow rate and shaft speed for each component. 
The compressor, turbine, and power turbine are assumed to be adiabatic, i.e., 
there is no heat transfer from these components to the environment. The 
turbocharger dynamics are controlled by the rotor inertia and dancing. The 
turbocharger turbine and compressor speeds are equal. Finally, the power 
turbine shaft is connected to the engine crankshaft via a specified gear ratio 
transmission. Hence, the power turbine speed can be determine- directly from 
the reclprocator engine speed. 
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3. CONSERVATION EQUATIONS 

In this section, equations for the conservation of mass and energy are 
developed for the contents of an open thermodynamic system. The conservation 
equation for the fuel mass is used to develop a differential equation for the 
change in fuel-air equivalence ratio of the system. The energy conservation 
equation is developed to obtain a differential equation for the change in 
temperature of the thermodynamic system. 

3.1 Conservation of Mass 

3.1.1 Conservation of total mass 

The rate of change of the total mass in any open system is equal to the sum 
of the mass flow rates into and out of the system: 


• • 

m - I m (3-1) 

J J 

Note tnat the convention used in our model assumes that mass flow rates 
into the open system are taken as positive, while mass flow rates out of the 
system are negative. 

3.1.2 Conservation of fuel ma ss 

In particular, conservation of the fuel species can be expressed as 


where ay denotes the fuel content in \e open system ( Includes fuel added by 
injection and fuel in the form of combustion products). 

Defining the fuel fraction, F, in the system as 

F » ay/m 

where m is the total mass in the system, equation (3-2) can be re-written as 
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d_ 

dt 


(mF) 




(3-3) 


where denotes the fuel fraction of the mass flow entering or leaving the 
open system through the j port. 

Differentiating the left-hand side of equation ( 3 — 3 ) gives 


raF 


Z m.F. - mF 
j J J 


( 3 * 4 ) 


or substituting for m from equation (3-1) results in a differential equation 
for the change in the fuel fraction of the open system, i.e. 


F - l (nij/m)(Fj- F) (3*5) 

An average fuel-air equivalence ratio, <t>, for the contents of the open 
system can be defined as 


Ul#*/ III 

0 - p^sf5 (3-6) 

Where m & is the mass of air in the open system and FASTO denotes the 
stoichiometric fuel to air ratio. 

Expressing the equivalence ratio, 0, in terms of the fuel fraction, F, 

i.e. 


* “ FASTO W (3 ’ 7) 

and differentiating (3~7) with respect to time we obtain an equation for the 
rate of change of the equivalence ratio of the open system, i.e. 


0 


__i F 

FASTO (1-F) 2 


(3-8) 


with F given from equation (3*5). 
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3.2 Conservation of Energy 

The general energy equation for an open thermodynamic system may be 
written as 

f • « • 

E - l m.h - Q - W (3*9) 

j J J ^ 

with the rate of change of the energy of the system being given by 

« ■ 3;<" h > - dt <pV> «- ,0) 

viiere 


[m.h. is the net rate of influx of enthalpy 
j J J 

• • 

Qy- Z is the total heat transfer to the walls, i.e. the sum of the 
1 

heat transfer rates to the different surfaces of the control 
volume of interest 
• • 

W - pV is the rate at which the system does work by boundary 
displacement. 

The dots denote differentiation with respect to time. Note that the 
convention used is that heat loss from the system and work done by the system 
are taken as positive. 

Differentiating the left-hand side of equation (3-9) gives 

• • 

mh - I m.h - 
j J J 

The contents of any open system, i.e. air and combustion products, can be 
represented as one continuous medium by defining an average equivalence ratio 
at each point in time. Gas properties are obtained assuming ideal gas 


V p V - 


mh 


(3-11) 
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behavior and thermodynamic equilibrium (see Appendix A). With these 
assumptions, we can express the enthalpy, h, and the density, p, of the 
mixture of air and combustion products as 

h - h(T,p,<|>) (3-12) 

P - p(T ,p, 4>) (3-13) 

Hence, the rate of change of the above fluid properties with respect to time, 
or crank-angle, can be written as 


Where 


h - c T ♦ c T p ♦ c 
P 1 <P 


„ r 311 ^ 

p" 3T ) p,* 


V 

c * (--) 

♦ Wt,p 


and 


p „ (§£) x + ( 3 £) p + ( 3 fi) 

p v arp,4. 3p ; T,4> p v 3<i> T,p 
The equation of state for ideal gases 
p - RpT 

can be expressed in differential form as 


(3-1*0 


(3-15) 


( 3 - 16 ) 



p R p T 


Re-arranging equation (3-17) and using (3-16), we can write 



(3-17) 


( 3 - 18 ) 


Substituting for p from equation (3-15) into the above equation, we can 
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express R in terms of p, T and $, i.e. 


r a (! B. 2fi)p - (-B- ♦ -2- ^£)x — 2- -2 a 

" V VP O WP ' 9 O JT' 1 O !»A ▼ 


pT p^r* V p^ 31 p^r 3 * 

From the differential derm of the equation of state, we can express the 
time rate of change of pressure as 


• • 


„ ,R m T V, 

p “ p *R m f”v^ 


(3*20) 


or substituting for R from equation (3—19) * and with sane manipulation, we 
Obtain the time rate of change of pressure: 


; . -IL- ( . V . 1 • . 1 fe • m, 

P 3p/3p V V p 3T par to 


(3-21) 


Returning to the energy equation (3~1 1 ) , and expressing h in terms of 
its partial derivatives with respect to T, p and <p, we obtain 


roc, 


t p T - £ nijhj- Qy+ (V - m^Jp - me^ <p 

J 


- hra 


(3-22) 


Finally, substituting for p from equation (3-21), we obtain an 
equation for the rato of change of temperature that does not explicitly depend 
on the rate of change of pressure of rthe system, i.e. 


me 


'j - « Vj- V _ ilS **•> (3 - 23) 


‘* M "> 8 * (575) 0 ' "V 


(3-2«) 




'H 


V 


i 
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Dividing (3-23) by m, and collecting terms in T, <|> and m, we get 

‘V ! to* ■ i ( j "j h J- V - <V ! U> ♦ * 6 5 <’ - §> - B l (3-25) 

leading to 


§r™/i h, v C ; 1 


T = A C m ' B } * V ' B ♦ + M iz m j h j“ V ] 

J 


where 


(3-26) 


A - c ♦ § U - o ♦ ^I2 ( l . c ) 
p p 3T p ( dp/dp) ( p V 

c - C + 5 & . c * i*£C*ihl . c J 

♦ P <J> Op/3p)'p C t' 


(3-27) 

(3-28) 


3.3 Application of Conservation Equations to Reolprocator Cycle 

The above conservation equations can be applied to the thermodynamic 
processes of the four-stroke diesel cycle as follows: 

1) Intake 


m 


m. - 
in 


m 

ex 


(3-29) 


F “ lT< F in- F > ' V< F ex- F > 


(3-30) 


• • 

t - B r m /i hv V C * 1 .* 

A C m (1 " B ) ' V * B ♦ + BiS (m in h in‘ m ex h ex' V 


(3-31 ) 


11) Compression 

(3-32) 
(3-33) 


m rn 0 
F - 0 
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(3-3U, 




i34s5 




T - -[ "] 

1 A L V Bro J 


(ill) Combustion 


m * id. 


( 3 - 35 ) 


• m f 

F - "(1 - F) 

ID 


( 3 - 36 ) 


BrV 


T • - s> - M ♦ ♦ k<Vf- 


( 3 - 37 ) 


Where h f is the absolute fuel enthalpy, given by (B-1 ) . 
(iv) Exhaust 


m 


m 


ex 


( 3 * 38 ) 


m 

- -”(F - F) 
m ex 


( 3 - 39 ) 


rn 


T B r ‘“ex,, h, V C * 1 ’ . * 

T " jf m~ ~ B ” V ~ B ^ Bm ™ex h ex~ V 3 


( 3 - 40 ) 


Note that the fuel fraction of the mass flow rate through the exhaust 


port, F , could be different from the cylinder fuel fraction, F, in reverse 

wv 

flow situations. 
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4. MODELING OF RECIPROCATOR ENGINE PROCESSES 
4.1 Gas Exchange 

A one-dimensional quasi-steady compressible flow model is used to 
calculate the mass flow rates through the intake and the exhaust valves during 
the gas exchange process. The intake manifold and the exhaust port are 
treated as plenums with known pressures. Furthermore, the temperature and 
average equivalence ratio of the intake charge (fresh air at intake manifold 
conditions) and the exhaust charge (mixture of air and combustion products at 
cylinder conditions) are known. When reverse flow into the intake manifold 
occurs, a rapid-mixing model is used, i.e., instantaneous mixing between the 
back-flowing charge and the intake charge is assumed. 

At each step of the gas exchange process, values for the valve open areas 
and discharge coefficients are obtained from tabulated data. Given the open 
area, the discharge coefficient, and the pressure ratio across a particular 
valve, the mass flow rate across that valve is calculated from: 


P D 2 P s 2/Y 

°d A Rf VYRf o { YH[ ( p') 
o o 


p (Y+1)/Y 1/2 

(p § > 

o 


where 


Cj « discharge coefficient 
A - valve open area 

p Q - stagnation pressure upstream of restriction 
p - static pressure at restriction 

3 

T q • stagnation temperature upstream of restriction 
Y » ratio of specific heats 
R - gas constant 


(4-1) 



When the kinetic energy in the cylinder is negligible, the stagnation pressure 
and temperature reduce to the static pressure and temperature, respectively. 
For the case of choked flow, equation (4-1 ) reduces to 

p (Y+1)/2(Y-1) 

• ■ Vn 7Sf o < w ) 

0 

The mass, m(t), in the cylinder at any time t can be found from 
integration of the mass conservation equation (3*1 ) , i.e. 

m(t) - V £ m in dt ' t m ex dt iu ~ 3) 

o o 

where m Q is the mass in the cylinder at the start of the cycle calculation, 
inlet valve opening (IVO). 

4.2 Combustion Model 

The diesel combustion process is a complex, unsteady, heterogeneous, 
three-dimensional process. A complete mathematical combustion analysis would 
require accurate models of compressible viscous air motion, fuel spray 
penetration, droplet break-up and evaporation, air entrainment into the spray, 
combustion kinetics, turbulent diffusion and so on. The details of the 
combustion process would depend on the characteristics of the fuel, the design 
of the combuotion chamber and the fuel injection system, and on the engine’s 
operating conditions. Although an adequate conceptual understanding of diesel 
combustion has been developed t . oate, a comprehensive quantitative model of 
all the individual processes has yet to be proposed. 

A useful approach to the problem of combustion simulation has been to 
model combustion as a heat release process, as originally proposed by Lyn [7]. 
The rate of heat release (or, equivalently, the rate of fuel burning) can be 
defined as the rate at Which the chemical energy of the fuel is released by 
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combustion. Based on heat release analysis and high-speed photography, Lyn 
has provided an excellent description of the different stages of diesel 
combustion. These stages can be identified on the typical rate of heat 
release diagram for a DI engine shown in Fig. 2 as follows: 

Ignition delay : The period between the start of injection (the point 
at Which the injector needle starts to lift) and ignition (the start 
of positive heat release due to combustion). 

Pre-mlxed or rapid combustion phase : In this phase, combustion of 

the fuel which has mixed with air to within the flamability limits 
during the ignition delay period occurs rapidly in a few crank angle 
degrees. This results in the high initial rate of burning generally 
observed in direct-injection diesel engines. 

Mixing controlled combustion phase : Once the fuel and air premixed 
during the ignition delay have been consumed, the heat release rate 
(or burning rate) is controlled by the rate at rfiich mixture becomes 
available for burning. The heat release rate may or may not reach a 
second (usually lower) peak In this phase; It decreases as this phase 
progresses. 

Late combustion phase : Heat release continues at a low rate well 
into the expansion stroke. Eventually, the burning rate asyptoti- 
cally approaches zero. The nature of combustion during this phase Is 
not well understood. Possible processes are that a small fraction of 
the fuel may not yet have burned, or energy present in soot and fuel- 
rich combustion products can still be released, etc. Given the some- 
what arbitrary limits of this phase, combustion models usually focus 
on the main heat release periods, l.e. the pre-mixed and mixing 
controlled combustion phases. 
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Based on his observations, Lyn [7] proposed a method of predicting 
burning rates from tne rate of fuel injection. The fuel injected is divided 
into elements according to the order in which they enter the combustion 
chamber. These fuel elements become "ready for burning" according to a 
certain law. Thus, a "ready for burning" diagram can be obtained from the 
rate of injection diagram. At the ignition point, which occurs after the 
lapse of the delay period, the part of the injected fuel which has already 
been made "ready for burning" is added onto the current "ready for burning" 
diagram, causing the initial Sharp peak in the burning rate diagram. Subse- 
quent burning is essentially governed by the rate of injection. Although this 
method gives a reasonable fit to data obtained over a wide range of speeds and 
loads, it requires further refinement and calibration before it can be used in 
computer simulation work. 

An alternative approach to modelling combustion, in the context of 
computer simulations predicting engine performance, is to describe the 
apparent fuel burning rate by algebraic expressions. The constants in these 
expressions can be chosen suitably to reflect the dependence of the actual 
fuel burning rate on engine type and particular operating conditions. 

Shipinski et al [8] attempted to correlate the apparent rate of fuel 
burning with the rate of fuel injection by fitting a Wiebe function to heat 
release diagrams obtained from tests on a high-speed swirl-type direct 
injection engine. Although Shipinski obtained a reasonable agreement with his 
engine data, the heat release shape defined by the Wiebe function alone has a 
notable difference from the two-part characteristic with the initial spike 
that is measured on most engine types. 

To overcome tills problem, Watson et al [9] proposed that the apparent 
fuel burning rate could be expressed as the sum of two components, one 
relating to pre-mixed and the other to diffusion-controlled burning, l.e. 


(4-4) 


BL 


m + 
P 


where m Is the fuel burning rate with respect to crank angle and subscripts 
t, p, d denote total, pre-mlxed and diffusion burning, respectively. 

In order to quantify the proportion of the fuel burnt by either 
mechanism, a phase proportionality factor, B, Is Introduced. This expresses 
the cumulative fuel burnt by pre-mlxed burning as a fraction of the total fuel 
injected, i.e. 


S 



(4-5) 


Consequently, a non-dimensional apparent fuel burning rate curve can be 
written as 


MJt) - BM (i) + (1 - B)M.(t) (4-6) 

t p d 

where M(t) is the non-dimensional burning rate distribution, and t the 
normalized crank position. 

The phase proportionality factor, B. Is considered to be controlled by 
the length of the ignition delay period (since the fuel that Is prepared for 
burning during this period governs the pre-mlxed burning phase), and the 
overall cylinder equivalence ratio, a . This can be expressed by a relation 
of the form 

♦ b 

6 - 1 - a-~- (H- 7 ) 

ID C 

there ID is the Ignition delay (see Section 4.3), and a, b, c are suitable 
constants to fit engine cylinder pressure data. 
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Fut:ienrore, Watson concluded that the best representation of the 
experimental data was achieved using the following ^.*)onent burning rate 
distributions: 

Premixed burning : 


M p (T> - 1 - (1 - t Cp1 ) Cp2 


(CL.-1) C nl (0,-1) 

OT M n (T) “ C n1 C n? T ? 0 " t P > P? 

P pi P2 


('•-8a ) 


(4-8b) 


Diffusion controlled burning (Wiebe function) 

M < j(t) - 1 - exp(-C d1 T j (4-9a) 

or M d (x) - C^C^x a ' exp(-C d1 t “-) (4-9b) 

where C p1 , C^, C d1 , are shape factors. The shape factors in equations 
(4-8) and (4-9) can be determined as a function of the engine operating 
conditions. 

Using data from three typical turbocharged truck engines, (design details 
are summarized in Table 1 ) , Watson established that the pre-mixed burning 
factors Cp^ and were adequately modeled, for all three designs , by the 
following correlations: 

C p1 - 2.0 ♦ 1.25 x 10‘ 8 (ID x N) 2 *** (4-10a) 

c p2 - 5000 (4-10b) 

where ID ■ ignition delay, ms (see Section 4.3) 

N - engine speed, RPM 
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Engine 

Engine 

Engine 

Source: 
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TABLE 1 

Appropriate Heat Release Profile Constants 
For Different Engine Designs 


ENGINE 

A 

B 

C 

a 

0.296 

0.95 

0.81 

b 

0.37 

0.41 

0.28 

c 

0.26 

0.28 

0.51 

k i 

14.2 

16.67 

7.54 

k 2 

0.644 

0.6 

0.65 

k 3 

0.79 

1.2 

0.93 

k 4 

0.25 

0.13 

0.22 


A: 6-cylinder, in-line, turbocharged, 4 stroke, DI, with 

inlet port swirl generation, and deep bowl combustion 
chamber. 

B: V8, turbocharged and intercooled, DI, A stroke, diesel 
engine. 

C: Turbocharged and intercooled, but more highly rated 
than engine A or B. 


Ref. [10] 
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The difusion burning factors C d1 , controlling the rate and duration of 
diffusion burning, and C d2> controlling the timing of the peak diffusion 
burning rate, were modeled by correlations of the form 


C ... - k. <f> _k 2 

dl 1 ove 


C., - k. C *4 
dd 3 dl 


(4-1 la) 
(4-1 1b) 


where 4> ove is the overall cylinder equivalence ratio, and k 1 , k^, k^, k^ are 
constants appropriate for each combustion chamber design. Table 1 summarizes 
the values for the diffusion-burning constants in equation (4-11) and the 
phase-proportonality constants in equation (4-7) that Watson etablished for 
the three engines used in his tests. 

In our engine simulation, we have followed Watson's approach to describe 
the heat release profile. As a starting point, Watson's constants for Engine 
A were used to describe our engine. However, at full load and rated speed 
conditions, the calculated peak cylinder pressure and the brake mean effective 
pressure were correspondingly found to be 26% and 6.6% higher than test data 
ootained at Cummins. In order to calibrate these constants for our engine 
design, an iterative procedure was followed [10], where the shape factors C dl 
and were modified until the calculated cylinder pressure diagram was in 
good agreement with the pressure data that we obtained from Cummins (one 
pressure trace, at rated load and speed, at a given injection timing). 

We concluded that the Engine A values for a, b, c, k 1 , k 2 and k^ were 
adequate, while would have to be close to 1.05. These values were used in 
all our subsequent calculations. 
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4.3 Ignition Delay Model 


The ignition delay period in a diesel engine was defined in section 4.2 
as the time (or crank angle) interval between the start of injection and the 
start of combustion. The start of injection is usually taken as the time when 
the injector neeole lifts off its seat (determined from a needle lift 
indicator). The start of combustion is more difficult to determine. It is 
best identified from the change in slope of the heat release rate versus time 
curve which occurs at ignition. 

Both physical and chemical processes must take place before the injected 
fliel can burn. The physical processes are: the atomization of the liquid 

fuel jet; the vaporization of the fuel droplets; the mixing of fuel vapor with 
air. The chemical processes are the precombustion reactions of the fuel, air, 
residual-gas mixture which lead to autoignition. These processes are affected 
by engine design and operating variables, and fuel characteristics. 

Ignition delay data from fundamental experiments in combustion bombs and 
flow reactors have usually been correlated by equations of the form 

ID - A p' n exp(E/RT) (4-12) 

where ID - ignition delay, ms 

E « apparent activation energy for the fuel autoignition process 
R - universal gas constant 
p - gas pressure, atm 
T - gas temperature, Kelvins 

and A and n are constants dependent on the fuel and the injection and air- 
flow characteristics. Representative values for A, n and E are given in Table 
2. Additional details can be found in [11]-[13]. 
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TABLE 2 


Constants for Arrhenius 1 
for Ignition Delay 


A n 



y a 






These correlations have usually been derived from tests In uniform air 
environments Where the pressure and temperature only changed due to the 
cooling effect of the fuel vaporization and fuel heating processes. However, 
in a diesel engine, pressure and temperature change considerably during the 
delay period due to the compression resulting from piston motion. Another 
problem is that the form of these simple correlations is not sufficiently 
flexible to allow all the influencing fuel and engine parameters to be 
included in the calculation of the ignition delay. 

Hardenfoerg and Hase [14] have developed an empirical formula for 
predicting the Ignition delay In DI engines as a function of fuel 
characteristics, engine parameters and ambient conditions. Dent [15] has 
shown that this formula can give reasonable agreement with experimental data 
over a wide range of engine conditions. However, the pressure and temperature 
used in this correlation are identified as the corresponding conditions at top 
deal centre, estimated using a polytropic model for the compression process. 

It is felt that such a poly tropic model is not appropriate in our simulation 
context, vhere pressures and temperatures can be accurately predicted over the 
duration of the ignition delay period. 

We have included two approaches for ignition delay in this simulation: 

a. The crank angle at start of combustion can be specified: this often 
is useful when it is desired to suppress changes in combustion timing 
(which may shift the oombustion process relative to its optimum location 
in the cycle). 

b. The ignition delay period is calculated as the difference between the 
time at Which combustion occurs (t^) and the time at Which injection 
starts (t in j) . The time t ign is Obtained by integrating over the 
duration of the ignition delay period the reciprocal of the predicted 
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ignition delay at each time step until the following relationship is 
satisfied: 


Agn dt 

' I6(t) * 

inj 


(4-13) 


where the instantaneous estimates of the ignition delay are calculated from an 
equation of the form of (4-12), with pressure and temperature taken at their 
instantaneous values. 


4.4 Heat Transfer 

The heat transfer mechanisms in a diesel engine include forced convection 

(Q c ) from the turbulent flow in the cylinder to the combustion chamber walls, 
and radiation (Q^) from the flame and the burning soot particles. The total 
heat transfer rate (C^) is therefore given by 

% - % * % ('•-I'O 

4.4,1 Convective heat transfer 

The convective heat transfer at the gas-to-cylinder wall interface will 
depend on the temperature gradient in the boundary layer at the surface. 
However, due to the inherent difficulties in calculating the details of 
turbulent fluid notion in the combustion chamber during the operating cycle of 
the engine, the convective heat transfer rate is usually expressed as 

Q c - h A (T g - T w ) (4-15) 

where 

h - convective heat transfer coefficient 
A » surface area 
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T • bulk mean gas temperature 
8 

T - inside wall surface temperature of cylinder head, piston, or 
w 

liner, as appropriate. 

»he problem is then to devise a method to calculate the convective heat 
transfer coefficient that appears in (4-15). The approach usually taken is to 
calc.: ’at 3 h from a Nusselt-Reynolds number correlation analogous to that used 
for steady turbulent flow in a pipe [l6]-[20], i.e., 

Nu = a Re d Pr e (4-16) 

where 

Nu - hL/A : Nusselt nunber 

Re - VL/v : Reynolds nunber 

Pr ■ we /A: Prandtl nunber 
P 

L - a characteristic length 
V - a characteristic velocity 
A - thermal conductivity 
v » kinematic viscosity 
u - dynamic viscosity 
p « density 

0 ^ - specific heat at constant pressure 
and a, d, e are constants adjusted to fit experimental data. 

Fortunately, there is little variation in the Prandtl number for air and 
combustion products, vbich is usually close to unity. Consequently, we may 
drop the Prandtl number dependence in equation (4-16) with little loss in 
accuracy, so that 

Mu - a Re d (4-17) 
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To calculate the convective heat transfer coefficient from correlation 
(4-17), instantaneous values for the characteristic length and velocity 
scales, and the gas transport properties (u, p, and X) are needed. Currently, 
it is not possible to predict these parameters and their spatial and temporal 
variation with arty accuracy in an internal combustion ergine. To overcome 
this difficulty, representative values of the characteristic length and 
velocity scales, and the gas temperature, pressure, and equivalence ratio at 
Which the gas properties are to be evaluated are chosen. 

For our heat transfer model, the characteristic ler^th scale is taken to 
be the macroscale of turbulence, as defined by equation (4—36) . The 
characteristic velocity V is postulated to be an effective velocity due co 
contributions from the mean kinetic energy, the turbulent kinetic energy and 
piston motion, i.e. 


V - [U 2 * u ,2 + (V /2) 2 ] 
P 


where 


(4-18) 


U - mean flow velocity, defined by (M-28) 
u' • turbulent intensity, defined by (4—29) 

V p - instantaneous piston speed 

While thiS expression for V is speculative, it is constructed in such a way 
that increases in any of the three component velocities lead to increases in 
the heat transfer rate, vdiile at the same time errors due to overestimating 
the contribution from any one component are minimized. 

Many attempts have been reported to determine the constants a and d, 
through curve-fitting experimental results [l6]-[20]. Suggested values are 


a - 0.035 to 0.13 
d - 0.7 to 0.8 


(4-19) 
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depending on intensity of change notion. The gas density and the transport 
properties, y and A, that appear in correlation (4-17) are evaluated at the 
mean gas temperature, pressure, and equivalence ratio of the cylinder contents 
(see Appendix C). 

4.4.2 Radiative heat transfer 

The primary sources of radiative heat transfer in a diesel engine are the 
high temperature burned gases and the soot particles which are formed as an 
intermediate step in the turbulent diffusion-controlled diesel flame. Because 
the particle size distribution, rumber density, and temperature, and flame 
geometry are not well defined in a diesel engine, direct measurements of 
radiation on operating engines are required. 

Estimates of the relative importance of radiation in cooled diesel 
engines have varied between a few and 50 percent of the total heat transfer 
[17],[l8],[21]-[28]. The limited experimental engine radiation measurements 
to date are summarized in Table 3* In general, the radiant heat flux depends 
on the location on the combustion chamber surface, crank angle in the 
operating cycle, engine load, engine size, and engine design. At high load, 
the measurements suggest that the radiant heat flux is betwseen 25 and 45 
percent of the total heat flux. 

Due to the complexity of the problem, accepted prediction formulas for 
radiant heat flux in a diesel engine are not available. Annand [17] has 
proposed a radiation term based on the average bulk gas temperature of the 
forms 


Or * - tJ 4 ) (4-20) 

where k » empirical radiation constant 
r 

A ■ surface area 
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TABLE 3 



Relative Importance of Radiant Heat Flux 


• » • # 




VL 

V Q w 


Engine 

Load Range 

Peak Values 

Mean Values 

Ref. 

DI, M stroke 

Mid-Full 

9-15 

10-30 

[21] 

DI, 2 stroke 

Full 


35 -J t5 

[23] 

Prechamber 

Idle-Full 


7-23 

[24] 

DI, swirl 

Light-Full 


0-40 

[25] 

DI, swirl 
4 stroke 

80*-Full 

12-18 


[26] 

DI, swirl 
4 stroke 

Light-Full 

70-13 

21-14 

[27] 

[28] 
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k - Co 
r r 


where o is the Stephan-Boltzmann constant (56.7x10 12 kW/m 2 ^) , and is 


T - average bulk gas temperature 
© 

T w - inside wall surface temperature of cylinder head, piston or liner, 
as appropriate. 

During the intake, compression and exhaust processes, the radiative heat flux 
is taken to be zero. During combustion, Annand and Ma [18] suggested that 

(4-21) 
an 

adjustable calibrating constant with values in the range of 0. 6-3.1, depending 
on the engine speed and load. Note that since the temperature used is the 
average bulk temperature and not the flame temperature, C, is not an 
emissivity. Limited evaluation of this approach has shown that » 0.6 gave 
approximately correct radiant flux magnitude for one engine but was too low 
for another [21], while C r - 1 .6 gave radiant heat fluxes higher than 
experimental data [22]. Flynn et al [28] have measured the instantaneous 
radiant heat transfer in a direct-injection diesel engine under a wide range 
of engine operating conditions. A monochromator was used to measure intensity 
of radiation at seven wavelengths. By assuming a monochromatic emissivity law 
consistent with observed radiation from ve*y small particles [29] and 
integrating over all wavelenghts, Flynn et al. obtained an apparent radiant 
temperature, an apparent grey-body emissivity, and a total radiant heat 
transfer rate at each crank angle during the cycle. The results of Flynn’s 
work indicate that the apparent radiant temperature is much higher than the 
average bulk gas temperature. In fact, during the time of peak heat release, 
the apparent radiant temperature was found to be very close to the flame 
temperature. Furthermore, during the period of maximum radiation, the 
apparent emissivity was 0.8 to 0.9, and it dropped almost linearly to zero by 
the end of the combustion process. 
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In the light of Flynn's results, and due to the absence of any real 
fundamental basis or experimental support for Annand's model, an alternative 
radiation model was developed for the present work. The instantaneous radiant 
heat flux Is expressed as 

Q„ - e oA(T 1 4 - T **) (4-22) 

r a r w 

v/here e - apparent grey-body emissivity 
a 

o * Stephan-Boltzmann constant 
A - surface area 

» apparent radiant temperature 

Flynn's data suggest that the apparent radiant temperature Is close to 
the adiabatic flame temperature during the period of peak heat release. The 
adiabatic flame temperature can be modeled as the temperature of slightly 
greater than stoichiometric zones of hydrocarbon-air combustion products, 
i.e., T($ » 1.1). However, as combustion progresses and relatively fewer 
close-to-stoichlometric fuel-air zones are found in the cylinder, this 
adiabatic flame temperature becomes considerably higher than Flynn’s apparent 
radiant temperature [28]. A better estimate (in reasonable agreement with the 
data) of the apparent radiant temperature was found to be the mean of the 
adiabatic flame temperature and the average bulk gas temperature, i.e. 

T ♦ T(* - 1.D 

T p - - S g (4-23) 

The temperature of combustion products at an equivalence ratio of 1.1, 

T(<b - 1.1), is computed as a function of the instantaneous air temperature and 
pressure from a correlation obtained by applying a least-squares curve-fitting 
technique to results of the NASA equilibrium program [30], for constant 
pressure hydrocarbon-air combustion. Satisfactory accuracy (less than It 
error) was obtained by considering two adjacent air temperature ranges, i.e.: 


3 ^ 





T($ - 1.1) - [1 ♦ 0.000231 7 (T air - 950)] x (2726.3 ♦ 0.9306p - 0.003233p 2 ) 
for 800 K < T alp < 1200 K (4-24) 

T(+ - 1.1) - [1 ♦ 0.000249 (T alr - 650)] x (2497.3 ♦ 4.7521p - 0.11065p 2 ♦ 

♦ 0.000898p 3 ) 

for 450 K < T , < 800 K. (4-25) 

a *r 


The instantaneous air temperature, T alr , is calculated assuming adiabatic 
compression of the air from the condition at the start of combustion 
(subscript ign), i.e. 

T air- T air,ign ( ^ign )(Y " 1)/Y (4 ' 26) 

where 1 is the ratio of specific heats for air at the instantaneous 
temperature and pressure. 

The apparent emissivity is assumed to vary linearly between its maximum 
value (taken as 0.9) and zero over the duration of the combustion process, 
i.e. 


e (t) - 0.9(1 
a 



where t is the time vhen the exhaust valve opens, 
evo 


(4-27) 


4.5 Turbulent Flow Model 

The heat transfer model of the cycle simulation requires estimates of the 
characteristic velocity and length scales. To estimate these scales in a way 
«h!oh incorporates the key physical mechanisms affecting charge motion in the 
cylinder, a turbulent flow model Is used. This model is a variation of the 
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models used by Mansouri et al [ 31 ], Poulos and Heywood [32] In previous engine 
simulation work. 

The turbulence model consists of a zero-dimensional energy cascade. Mean 
flow kinetic energy, K, is supplied to the cylinder through the valves. Mean 
kinetic energy, K is converted to turbulent kinetic energy, k, through a 
turbulent dissipation process. Turbulent kinetic energy is converted to heat 
through viscous dissipation. When mass flows out of the cylinder, it carries 
with it both mean and turbulent kinetic energy. Figure 3 illustrates the 
energy cascade model. 

At any time during the cycle, the mean flow velocity, U, and the 
turbulent intensity, u* , are found from knowledge of the mean and turbulent 
kinetic energies, K and k, respectively. Thus, the following equations apply: 

K - ^ m U 2 (4-28) 

k • | m u ’ 2 (4-29) 

where the factor 3 in equation. (4-29) comes from assuming that the small scale 
turbulence is isotropic (and accounting for all three orthogonal fluctuating 
velocity components). 

The time rate of change of the mean kinetic energy, K is given by 


dK I' 2 . „ "e 

dt - 2 “i V i ' P " K u 


(4-30) 


Similarly, the rate of change of the turbulent kinetic energy, k, is 


* 

dt 


• P - me 



with 


c „ «’ 3 . LOOti l ^ 2 

c 1 l 


(4-31 ) 
(4-32) 
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where in » mass in the cylinder 


m,« mass flow rate into the cylinder 

A. 


m e » mass flow rate out of the cylinder 

Vj- jet velocity into the cylinder 

P - rate of turbulent kinetic energy production 

e - rate of turbulent kinetic energy dissipation per unit mass 

A * rate of turbulent kinetic energy amplification due to rapid 

distortion 

l « characteristic size of large-scale eddies 
In equations ( 30) and (4-31 ) , the production term P has to be defined in 
terms of flow and geometrical parameters of the chamber. However, since the 
above model does not predict spatially resolved flow parameters, P must be 
estimated from mean flow quantities only. 

Assuming that turbulence production in the engine cylinder is similar to 
turbulence production in a boundary layer over a flat plate [33], we can 
express P as 


(«) 

"tV 


P ■ M£> 


(4-33) 


2 ,. 


wh*»re p .■ C k /(me) is turbulent viscosity, and C - 0.09 is a universal 
t u u 

co istant. Again, as the velocity field In the cylinder is not known, the 

velocity gradient (3U/3y) is approximated as 


(^) 2 - c fr* 

l 3y ; Vl/ 


(4-34) 


where C. is an adjustable constant and L is a geometric length scale. 
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Using equat: ; ons (M— 33) , 34) , (4-28) and (4~32), we can express P as 


o 3/2 ki k 1/2 

2(|) C Cj*f)Q 
2 u 0 ^2 m 


(4-35) 


Furthermore, the characteristic size of the large-scale eddies, l, and the 
representative geometric length scale, L will be both identified with the 
macroscale of turbulence, assumed to be given by 


l - L - V/(vB 2 /4) (4-36a) 

where V is the instantaneous volume of the combustion chamber and B is the 
cylinder bore, subject to the restriction that 

L £ B/2 (4-36b) 

Hence, equation (4-35) can be re-written as 

K k ^ 2 

P = 0.3307 C g (*)(*) (4-37) 

During the compression and the combustion processes, the turbulent 
kit., tic energy decays due to viscous dissipation. At the same time the 
turbulent kinetic energy is anplified due to the rapid distortion that the 
cylinder charge undergoes with rising cylinder pressures. Consequently, an 
amplification term. A, was added to equation (4-31) to account for this 
effect. The anplification term will be larger during combustion vfcen the 
unburned gas is assumed to be compressed by the flame at a sufficiently high 
rate. However, in the diesel engine context of high compress ion ratios, the 
anplification term is included during compression, too. 

Using equation (4-29), the rate of turbulence amplification due to rapid 
distortion can be expressed as 


A 


du* 

3 m u dt 


(4-38) 
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where the rate of change of the turbulent intensity du'/dt can be estimated 
assuming that conservation of mass and angular momentum can be applied to the 
large seal eddies during the rapid distortion period. 

Under these assumptions, conservation of mass for a single eddy of volume 
V. requires that 

U 


p V , - p V, 

L o Lo 


(4-39) 


Where p is the mean gas density and subscript o refers to the conditions at 
the start of compression. Then, since 


V 


L 


(4-40) 


where L is the macroscale of turbulence, we can re-write (4~39) as 



Conservation of eddy angular momentum requires that 

U L - U L (4-421 

to wo o 

where U^ is the characteristic velocity due to eddy vorticity. 

Combining (4-41) and (4-42) with the assumption that 

U - u' (4-43) 

ti) 

the following relation is obtained for the evolution of the turbulent 
intensity during the rapid compression period 


.y: . (fi _) 

u* " 


1/3 


(4-44) 


o o 

Differentiating both sides of equation (4-44) and re-arranging we get 


Ho 


du* u’ dg 
dt * 3p dt 


(4-45) 



Hence, combining equations (*l— 38) and (4— *15) » the rate of turbulence 
anplification is given by 

A - ®- 2 Jf (4-46) 

p Qt 

or, introducing (4-29) • 

A - | k fi (4-47) 

3 P 

with p given from equation (3~15). 

4.6 Engine Friction Model 

To convert Indicated engine performance quantities to brake engine 
performance quantities, engine friction estimates are required. However, the 
measurement and analysis of engine frictional losses are yet to be 
satisfactorily resolved. This is primarily due to the inherent problem of 
direct, accurate measurement of these losses under actual running conditions. 
This problem occurs because the total loss is a summation of losses arising 
from the operation of the many components of the engine, and these components 
respond differently to changes in pressure, temperature and speed. 

Direct motoring of an engine is the common method of measuring losses, 
but clearly the motoring losses are not the same as the losses under firing 
conditions. Some of the reasons are the lower pressure acting on piston rings 
and bearings, the lower temperatures of the piston and cylinder bore surfaces 
and thus the greater oil viscosity, the greater running clearance of the 
piston, and the missing exhaust blowdown period. 

Nevertheless, a breakdown analysis of motoring losses supplemented by 
experiments on piston and ring friction rigs can be used to identify the 



relative importance of the many components of the total friction and their 
response to changes in design variables. In general, the rcxnponents of the 
losses expressed in terms of mean effective pressure, mep, tend to fall into 
three groups: 

(i) Losses due to boundary lubrication, where the friction forces are 
approximately invarient with speed. These losses are undoubtedly 
influenced by the compression ratio. 

(ii) Losses associated with hydrodynamieally lubricated surfaces in 

relative motion, which vary directly with speed. All major rotating 
parts fall into this group. 

(iii) Losses associated with fluid (air, water and oil) punping, which 
vary as the square of the speed. 

Therefore, the motoring losses can be expressed in the form 

F - A ♦ BN + CN 2 (4-48) 

vftere F are the losses in mep, N is engine speed, and A, B and C are 
constants. 

Millington and Hartles [34] have measured motoring losses on a large 
variety of automotive diesel engines during the course of development of 
prototype engines. Their work suggests a readjustment of equation (4-48) 
copied with suitable selection of the constants as follows: 


F 


* 7 -° Tooo 


♦ 1.5 


V 

(— — ) 
M000 ; 


(4-49) 


where F - motoring mep, psl 

A » compression ratio minus 4 for a DI diesel 
N - engine speed, RPM 
V - mean piston speed, ft/min 


42 


Zf 


. V 



Equation (4-H9) represents a sound empirical correlation of the motoring loss 
data obtained from diesel engines. It is used to obtain brake quantities from 
the indicated quantities computed in the engine simulation. 



5. MODELING OF OTHER SYSTEM COMPONENTS 


5.1 Turbomachinery Modeling 

Steady-state performance maps give the interrelationships among mass flow 

rate, efficiency, pressure ratio and rotor speed for each of the three 

turbomachinery components: turbocharger compressor, turbocharger turbine and 

power turbine. The map variables of mass flow rate, m , and rotor speed, 

map 

“map* are correc '' ec * by factors relating actual inlet conditions to standard 
conditions. The speed correction factor involves inlet temperature, and the 
mass flow rate correction factor involves inlet temperature and pressure, so 
that 


u - u> „ , (T . . /T ) 
map actual std in 


1/2 


(5-1) 


’ * 1 /? 

m -m . ,(p Vp. )(T. /T ) (5-2) 

map actual std in in std 

where the subscripts in and std refer to actual inlet and map standard inlet 
conditions, respectively. 

The turbomachinery maps, usually obtained in graphical form, are entered 
into the simulation in tabular form, i.e. as a one-dimensional array of rotor 
speeds, ranked in ascending order; and a three-dimensional array of the 
remaining map variables, arranged as follows: for each of the speeds stored in 
the speed array, a uniform number of data points is recorded, each consisting 
of the values of mas3 flow, efficiency and pressure ratio at that point. The 
data points for each speed curve for each map are ranked in ascending order of 
pressure ratio. Note that for the compressor map this implies entering data 
in descending order of mass flow. 

At a particular step in the cycle simulation, the tables need to be 
interpolated to find the necessary information. Appropriate routines were 




developed to interpolate the performance maps of the various turbomachine: y 
components. In general, these routines perform two-dimensional interpolation 
to calculate two unknown map variables from two known variables. Speed is 
always one of the two known variables and efficiency is one of the two unknown 
variables. Then, either mass flow rate is known and pressure ratio unknown, 
or vice-versa, depending on the turbomachinery component and the system 
configuration. 

The method of map interpolation, as applied tor example to the compressor 
map, is the following. Given a corrected speed, a search of the compressor 
speed array is performed, from the lowest value (i.e. the first array element) 
until a speed greater than the input speed is found. Using that greater speed 
and the speed just previous to that (the lesser speed), a speed interpolation 
parameter is calculated. This speed interpolation parameter is used to 
calculate values of pressure ratio at the input speed, until a pressure ratio 
greater than the input pressure ratio is found. This pressure ratio search, 
from low to high, will thus define a pressure ratio interpolation parameter. 
Then, using the speed and pressure ratio interpolation parameter, a two- 
dimensional linear interpolation is performed to calculate the mass flow and 
efficiency thit correspond to the input values of speed and pressure ratio. 

Certain important physical constraints mu3t be considered in modeling 
turbomachinery component performance: 

(a) Compressor Surge Line : When the mass flow rate through a compressor 

is reduced while maintaining a constant pressure ratio, a point arises at 
which local flow reversal occurs in the boundary layers. This will 
relieve the adverse pressure gradient until a new flow regime at a lower 
pressure ratio is established. Then, the flow will build up to the 
initial conditions, and thus flow instability will continue at a fixed 



frequency. Thi3 phenomenon is called surge. Clearly, a compressor 
should not operate in the low-efficiency, unstable region, to the left of 
the surge line. 

(b) Turbine Choking Characteristics : The mass flow range of a radial 

turbine, 3uch as the turbocharger and the turbocompounded turbines, is 
limited by choking at high pressure ratios. The choking characteristics 
of the turbine are speed dependent. This effect is caused due to the 
centrifugal field created by the speed of the rotor. 

To avoid any potential problems associated with the above constraints, certain 
provisions have been incorporated in the logic of the turbomachinery 
interpolation routines: 

(i) Since the pressure ratio versus turbocharger speed line i3 fairly 
flat close to the compressor surge limit, small changes in pressure ratio 
can result in disproportionately large changes in output mass flow rate 
values. In order to avoid any oscillations in the system convergence 
procedure vhich could result in a substantial increase in computational 
time, the pressure ratio versus mass flow rate curves are modeled as 
single-valued, with a small, non-zero, slope close to the surge line. 

(ii) During the turbocharger matching calculation, the rotor speed at 
some instant could become too low for the required boost pressure ratio. 
This means that the mass flow would have to be to the left of the ..urge 
boundary. Under these circumstances, the speed of the rotor is increased 
until an acceptable mass flow solution at the given input pressure ratio 
is obtained. 

(iii) During the course of the simulation iterations, certain engine 
operating conditions could correspond to points that lie beyond the 
normal operating regime of the turbomachinery maps. Then, linear 
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extrapolation la performed to extend the range of the map 
characteristics. This extrapolation is subject to certain checks, so 
that the map regimes are not extended beyond the turbine choking 
characteristics or to the left of the compressor surge line. 

Furthermore, if the final converged solution for the engine- turbocharger 
matching is well beyond the normal map regimes, the turbomachinery used 
is not appropriate for the given engine design and operating conditions. 
The calculation should be repeated with more suitable turbomachinery 
components (different machine sizes, or higher component efficiencies). 


5 . 2 Tu rbocharger Dynamics 

The rate of change of the mechanical energy of the turbocharger rotor, 
E t/c> depends on the difference between the power required to drive the 
compressor (negative), ai.a the power delivered by the turbocharger turbine 
(positive): 

* • # 

E w -W ♦ W„ .. (5-3) 

t/c compressor turbine 

The change in mechanical energy relates to the change in rotor speed according 
to the turbocharger dynamics, i.e., 

* ' 2 

E^ • Juu * Bu (5 _J 0 

where 

J - rotational inertia of turbocharger 
B - rotational damping of turbocharger 
uj - angular velocity 


1*7 


The compressor and the turbine powers can be expressed as: 


W - m (h, - h,) 

compressor cl 2 


(5-5) 


W. . . - m. (h 0 - h n ) 

turbine t 8 9 


(5-6) 


Solving for the change in speed <j gives: 


((m c (h 1 - h 2 ) + m t ( h g - hg) - Bu> )/Ju 


(5-7) 


Note that the enthalpy changes across the compressor and turbine can be 
calculated assuming compressible flow across the two turbomachines, i.e., 


u R/c 

b p 2 P 

(h r h 2 ) - ---c(p-) "13 

c 1 


(5-8) 


(h Q - h 9 ) 


R/c 

p 8, P 


^ C( 5:> “ 1] 


(5-9) 


5 . 3 Turbocharger Matching Procedure 

The simulation code is set up to analyze two different system 
configurations: i) a single-stage turbocharged diesel system and, ii) a 

turbocharged turbocompounded diesel system. Suitable turbocharger matching 
procedures have been developed for each case. These are summarized below. 

5.3.1 Single-stage turbocharged case 
The system configuration for this case is shown in Fig. 4. At a given 
instant, the values of the variables describing the state of the various 
system components are known (from the integration of the system governing 
equations over the previous time step). These include the intake and exhaust 
manifold pressures and the turbocharger rotor speed. Additionally, the 
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compressor inlet pressure and the turbine exhaust pressure are fixed, i.e., 
atmospheric pressure les3 intake air filter pressure drop and atmospheric 
pressure plus muffler pressure drop, respectively. By relating the compressor 
discharge pressure ;o the intake manifold pressure and the turbine inlet 
pressure to the exhaust manifold pressure, through suitable pressure drops, 
the pressure ratio across each turbomachine is determined. Hence, the 
compressor and turbine maps can be entered using the calculated pressure 
ratios and the rotor speed (same for both turbomachines) as input3. The 
output of the map interpolation routines determines the mass flow rate and 
efficiency for each component for the next time step. From these, the power 
required to drive the compressor , and the power delivered by the turbine are 
determined. Any excess power will result in a change in the rotor soeed 
according to tr.e turbocharger dynamics, i.e., Eq. (5-4). Finally, the values 
of the other state variables for the next time step will be determined from 
solution of the mass and energy conservation equations, where the compressor 
and turbine mass fluxes are taken from the output of the turboroach inery 
interpolation routines. 

5.3.2 Turbocharged turbocompounded case 
This case is relatively more conpllcated. For this system configuration, 
shown in Fig. 1, the Intake side, and hence the compressor map treatment is 
identical to the single-stage turbocharged case. However, on the exhaust 
side, althaigh the available pressure ratio can be defined as for the 
turbocharged case, the division of the pressure ratio between the two turbines 
is not known a priori. This calculation requires a special iterative 
procedure, based on continuity of mass flow through the exhaust system. At 
each time step, the calculation is started by assuming a mass flow going 
through the first turbine (taken as the value at the previous time step). 
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Then, the turbocharger turbine interpolation routine is entered with the mass 
flow rate and rotor speed as inputs. From the turbine pressure ratio returned 
by the routine, and the known exhaust manifold pressure, the first turbine 
exit pressur' is determined. The latter, and the system exit pressure 
(atmospheric plus muffler pressure drop), specify the pressure ratio across 
the power turbine. Then, the power turbine map can be interpolated with the 
pressure ratio and shaft speed (gear ratio times reciprocator speed) as input, 
thus defining the mass flow through the power turbine. An updated estimate of 
the mass flow through the first turbine can be then calculated based on the 
previous estimate and the power turbine mass flow. The scheme used involves 
under-relaxation to increase the stability and efficiency of the matching 
procedure. For a converged solution, the mass flow through the turbocharger 
turbine (plus flow through the wastegate, if appropriate) must equal the flow 
through the power turbine. Once the turbocharger turbine flow and efficiency 
are established, the system state variables for the next time step can be 
determined following the same procedure as for the single-stage turbocharged 
case. 






5. 1 * Intercooler Model 

The 1 tercooler, situated between the compressor discharge and the Intake 
manifold, serves to increase the density of the charge air by lowering its 
temperature. The intercooler is modelled as a heat exchanger of fixed area, 
over-all heat-tranfer ''^'■'‘lclent and cooling flow rate. The change in charge 
air temperature is de .- ,% ,.ned from the non-dimensional heat exchanger 
effectiveness, €, 


e 


(T, 


hi 


T h2 ) 


( V V 


(5-10) 
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where the subscripts h and c refer to the charge air flow (hot) and the 
coolant flow (cold) respectively; 1 and ? refer to inlet and exit conditions; 
and 

* compressor discharge temperature 
T^ * intercooler discharge temperature 
T « coolant inlet temperature (assumed to be fixed). 

C I 

The heat exchanger effectiveness is either known as a design parameter, 

or can be derived from graphical correlations that are available for the 

various typical heat exchanger configurations [35]. From the latter, e can be 

determined as a function of the capacity rate ratio, C . /C , and the number 

min max 

of heat transfer units, N. , where C . and C are respectively the smaller 

A min max 

and larger of the products of charge air and coolant flow rates with their 
respective specific heats, and 

"a • Au/C nin 

where A - heat exchange surface area (fixed) 

U ■ over-all heat transfer coefficient based on A. 

Assuming that C is much larger than C . , the expression for effectiveness 
max min 

reduces to the following simple form: 

e - 1 - exp(-N ft ) (5-11) 

5.5 Exhaust Manifold Model 

In most engine simulation programs of this type, the exhaust manifold is 
treated as a global open system [6]. A simple plenum is used with the exhaust 
gases from each cylinder flowing into the system and mixing with the all of 
the gases in the manifold. This approach has three primary disadvantages. 
First, the exhaust gases from a given cylinder are diluted by the gases in the 
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manifold in a manner that is not representative of a real manifold, where 
individual runners keep the gases from different cylinders separated during 
much of the flow path between the exhaust valve and the turbine. This 
dilution decreases substantially the gas temperature oscillations that occur 
in the manifold and at the turbine inlet. 

The second disadvantage of a common plenum model is that it is not 
possible to follow the change in gas temperature as the gases flow through the 
exhaust manifold, mix with the gases in each section, and lose heat to the 
surroundings. With a common control volume, only one temperature is used to 
represent the temperature of the gases in the exhaust manifold, ftgain, this 
is not a good representation of a real exhaust manifold where the gas 
temperature varies along the length of the path between the exhaust valve and 
the turbine inlet. 

The third disadvantage with a single control volume is that it is 
impossible to match the approximate volume, the inside surface area, and the 
cross-sectional area of a real manifold with only one set of dimensions for 
the model. The approximate volume is important for the accurate determination 
for the exhaust manifold pressure which in turn is an important factor in 
determining the turbocharger performance. The cross-sectional area and inside 
surface area are important for determining the temperature of the gases at the 
inlet to the turbocharger. With the single control volume approach, one is 
forced to compromise one or more of these geometric factors when specifying 
the dimensions of the manifold to be used for the simulation. 

To avoid these disadvantages with a single plenum model, the exhaust 
manifold is divided into a series of connected open systems. The manifold is 
considered to be a composite of ports, runners, and a common plenum section. 
The port section is taken to be the volume contained within the head of the 
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engine. The runner represents the section of the manifold outside the head 
Where the gases from one cylinder flow without mixing with the gases from any 
other cylinder. The plenum is defined a3 the volume where the gases from each 
of the engine cylinders, ports and runners are mixed before entering the 
turbine. A fourth system is also followed during the simulation that is 
composed of the runner from each cylinder and the plenum. This system has the 
properties that represent the average properties of the separate 3ub-systems, 
and it is used to determine the change in pressure in the manifold based on 
the overall mass balance. 

As is done with the engine cylinders, the simulation only calculates the 
properties of the ports and runners for one master port and one representative 
master runner. The master port and runner information necessary to include 
the other cylinders in the engine calculation is stored in high-speed memory 
and is retrieved with the appropriate phase shift to determine what was 
occurring in other ports and runners at any given time during the cycle. 


5.6 Manifol d C onser v ation Equation s 

The general manifold open system model is shown in Fig. 5. The 
differential equations for the change in total mass and fuel fraction for the 
manifold come directly from the continuity relations derived earlier in 
Section 3.1 , i.e., 


m 



( 5 - 12 ) 



(5-13) 
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where the subscript m refers to the intake manifold, or exhaust manifold 
section, and the subscript j refers to the j-th mass flow that enters or 
leaves the system of interest. 

The mass flows entering and leaving the intake manifold include the 
compressor mass flow (entering) Which is found from the compressor map; and 
the engine intake mass flow (leaving or entering) which is found from the 
reciprocator model. The mass flows entering and leaving the exhaust manifold 
as a Whole include the turbine mass flow (leaving) which is found from the 
turbocharger turbine map; the engine exhaust mass flow (altering or leaving) 
which is found from the reciprocator model; and the wastegate mass flow 
(leaving). 

The differential equation for the change in intake manifold temperature 
i3 derived from the generalized temperature equation applied to an open 
system, i.e., 

• 

T - k<*> (B - h m ) - ce + 1 (£ m.h. - Q)] (5-14) 

cn A in m m ta . j j 

m m j 

vhere A, B, and C are defined in Eqs. (3-27), ( 3~24 1 and (3-28), and are 

calculated by the thermodynamic property routines; Q is heat transfer rate to 

• • 

the manifold walls; and $ is related to F and F by Eq. (3-8). 

The differential equation for the change in intake manifold pressure is 
derived from the generalized pressure equation applied to an open system, i.e: 


P„ 


— e~ r- 1 t * (5) I 1 1 
3p/3p L p 3T m V m p 3<t> V 


(5-15) 


where the density derivatives are calculated by the thermodynamic property 
routines. 
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The exhaust manifold is broken down into a series of connected open 
systems that are contained within the overall manifold system (see section 
5.5). For these sub-systems, the pressure derivative is determined by applying 
the conservation equations (5-12) to (5-15) to the exhaust manifold as a 
whole. The exhaust ports are considered outside of the overall exhaust 
manifold control volume. However, the pressure derivative for the exhaust 
ports is set equal to that for the exhaust manifold as a whole. 

The rate of change of mass within each of these sub-systems can be 
obtained by rearranging equation (5-15) as follows: 


m = V (-^ p + -fi T + a ) 

i lap p i 3T l 3<J> V 


( 5 - 16 ) 


where the subscript i refers to the different sub-systems of the exhaust 
manifold. The rate of change of temperature and composition of each section 
are calculated from equations (5-1 4) and (5-13), respectively. 

The mass flow between sections of the exhaust manifold is determined by 
the mass flow through the exhaust valve (see section ^.1) and the rate of ma3s 
storage within any sections upstream of the section in question due to changes 
in the properties of those sections. For instance, the mass flow out of the 
port section and into the runner section is found by subtracting the change in 
mass in the port as found by equation (5-16) from the mass flow through the 
exhaust valve. 


5.7 Manifold H eat Tr ansfer 

With the incorporation of an intercooler, heat transfer from the intake 
manifold to the environment becomes small enough to be neglected. Heat 
transfer from the exhaust manifold, however, is significant. Heat transfer 
from the gas in the exhaust manifold to the environment involves a combination 
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of forced convective heat transfer from the gas to the inside walls, 
conduction from the inside walls to the outside walls and to the water Jacket, 
and natural (and probably forced) convection from the outer surfaces to the 
environment. 

To determine the heat transfer in the exhaust manifold, each section is 
considered separately. Within each section, the gas temperature, the heat 
transfer coefficient* and the inside wall surface temperature are assumed to 
be uniform. Both the gas temperature and heat transfer coefficient are 
allowed to vary with time while the inside wall surface temperature is taken 
as constant with time. The inside wall surface temperature can be either 
specified as an input, or calculated from a specification of the conponent 
wall structure, in the manna' described in sections 6.2.2 and 6.2.3. 

5.7.1 Port heat transfer 

The heat transfer in the exhaust port is highly unsteady. As the exhaust 
valve first comes open, a high velocity jet of high temperature gases sets up 
recirculation zones in the port [ 36 ] that result in a high heat transfer 
coefficient. When the exhaust valve is fully open, the flow resembles 
turbulent pipe flow and the heat transfer coefficient is diminished. Then, as 
the exhaust valve is going closed, there is another period tfien a narrow jet 
of gases enhances the heat transfer by setting up a recirculation zone. The 
valve open period is followed by a much longer period with very low mass flow 
and a correspondingly low heat transfer coefficient. 

In order to quantify the heat transfer in the exhaust port, the results 
of Caton [37] are applied. Based on fine wire temperature measurements of the 
port gas temperature in a spark ignition engine, Caton arrived at the 
following correlations for the heat transfer during different phases of the 
exhaust valve opening: 
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Valve opening phase ( 1/D < 0.2 ): 


Nu - 0.4 Re. 


(5-17) 


Valve open phase ( 1/D > 0.2 ): 


Nu = 0.023 C R C EE Re°' 8 Pr°' 14 


(5-18) 


Valve closing phase ( 1/D < 0.2): 


Nu « 0.5 Re. 


(5-19) 


Valve closed: 


Nu * 0.023 Re°‘ 8 Pr°' 1< 


( 5 - 20 ) 


where Nu = Nusselt number 


Pr = Prandtl number 


Re^ = Jet Reynolds number 


Re » Pipe Reynolds number, based on cycle-averaged mass 

flow 

V = Velocity of flow through exhaust valve 

J 

V - Pipe flow velocity 

i - Valve lift 


D - Valve diameter 


v ■ Dynamic viscosity for flow 


C D » Correction factor for surface roughness 
R 


C_ - Correction factor for entrance effects 
EE 


Some modifications were made to these results for the present study. A 
correction factor was added to the correlation for the valve open phase to 


account for the sharp bend for an exhaust port. This factor is the correction 
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for bent pipes discussed in section 5.7.2. Also, during the period when the 
valve is closed, the mass flow is based on the instantaneous mass flow between 
sections instead of on the average flow over the complete cycle, as was done 
by Caton. For a turbocharged engine, thi3 flow is not zero, in general, 
because gas flow is induced in inactive sections of the exhaust manifold by 
the pressure pulsations produced when other cylinders exhaust. 

5.7.2 Exhaust manifold 

Turbulent pipe flow correlations are usually applied to the exhaust 
manifold heat transfer. For fully developed pipe flow in a straight pipe with 
large temperature gradients, the heat transfer coefficient can be derived from 
the following experimental correlation that relates Nusselt, Reynolds and 
Prandtl numbers [35]: 

Nu = 0.023 Re°' 8 Pr°‘ 3 (5-21) 

In the exhaust manifold of an engine, not all of the qualifying conditions for 

this equation apply. Most significantly, the flow is not fully developed, due 

to the short length of the pipe, and the pipe is not straight. To correct for 

these differences, two correction factors are introduced, based on empirical 

results of studies of the effects of entrance length and pipe bends on pipes. 

The correction factor for entrance effects is based on the work of 
Boelter, Young, and Iversen [38]. They conducted experimental work on the 
variation of the heat transfer coefficient along the length of a tube starting 
at the inlet of the tube. Fitting a curve to the experimental results for a 
tube with an elbow at the entrance, yields the following equation for the 
local heat transfer coefficient: 


Nu x 


- 0.3 


Se • s: - 2 - 2 <5> 


(5-22) 
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valid for 3 < x/D < 14 

where x » entrance length 

Nu * local Nusselt number 
x 

Nu^ ■ Nusselt number for fully developed flow. 


In order to find the average correction factor for a section, equation (5-22) 
is integrated over the section and divided by the length of the section, 
leading to 


C EE “ 3J ° 


(L (L 

0.3 (L 2 ; ( V 

‘ 4“ L"* 


(5-23) 


i^iere D « diameter of the section 

L, - distance f-'orn the exhaust valve to the inlet of the section 


» distance from the exhaust valve to the 


The correction factor for bent pipe is based 
McLaughlin [39] and Rogers and Mayhew [ JO]. They 


correction factor for bent pipe is: 

0.1 

0.05 


Nu, 


"'BP Nu 


BP 


( R> Re 


outlet of the section 

on the work of Seban and 
found that the appropriate 


(5-24) 


where r 
R 



Pipe radius 
Bend radius 

Nusselt number for the bent pipe 
Nusselt number for a straight pipe. 


The Nusselt nurriber for the bent pipe is the average of the heat transfer 
coefficient around the circumference of the pipe. 

Combining the formula for fully developed straight pipe flow with the 
correction factors for entrance effects and for pipe bends gives the 


61 




% U L ■ U lb ' *■ ^ ■ . -■ ' . 

i 

.' S 

/ correlation used for the heat transfer coefficient for sections of the pipe 

V- downstream of the port section, and for the exhaust port during the periods 

when the exhaust valve is fully open or closed: 

Nu - 0.023 C EE ,C Bp Re°* 8 Pr°* 3 (5-25) 

V 

V. 

.j 5.7.3 Connecting pipe between the turbines 

^ The heat transfer for the connecting pipe between the turbocharger and 

v the power turbines is treated in the same way as for the exhaust manifold. 

The heat transfer coefficient is based on turbulent pipe flow using equation 
(5-25). The overall heat transfer coefficient and the insiie wall surface 

C* 

M 

! temperature are calculated in the manner described in sec <ns 6.2.2 and 

•i 6.2.3* 

'a, 

I 

•1 

5.8 Pressure Losses 
' 5.8.1 General 

Pressure loss terms have been included at five locations in the overall 
system model. These are: 

between compressor discharge and intercooler inlet, 
across intercooler, 

between exhaust manifold and turbine inlet, 

* 

between turbine outlet and power lurDine inlet, 
between power turbine outlet and atmosphere. 

Each of these pressure drops is calculated using the corresponding friction 
factors and friction coefficients for the geometry of each passage. For 
straight-sections: 

AP - 4f(L/D)(pV 2 /2) (5-26) 

Where l, - length of passage 
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D « d latter of passage 
p » bulk density 
V - average velocity 

f - friction factor, which is correlated by (5-17) for the surface 
roughness and range of Reynolds numbers to be encountered. 

Note that 

f - 0.046/(Re)°' 2 (5-27) 

with Re » pVD/y 

For bends, enlargements, contractions, etc: 

AP - K f pV 2 /2 (5-28) 

where K^. - friction coefficient for a particular passage geometry. Values 
of Kj, for typical geometries are commonly available [35]. 

5.8.2 Exhaust man if old and turb in e connecting pipe 

The pressure drop for the exhaust manifold of a 6 cylinder diesel engine 
was found by Primus [41] to be a factor of 10 to 15 times higher than the 
pressure drop calculated based on equations (5-26) and (5-27) alone. This 
increase in flow losses is apparently due to the complex shape of the rr.nifold 
and the interaction of the flow with the open passageways from other 
cylinders. For this reason, equation (5-28) is used for the calculation of 
the exhaust manifold pressure drop with K^, left as an input parameter for the 
user to specify. Primus found that values of 2 to 3*5 were appropriate for 
his test manifold. Different values may be used for manifolds of different 
designs. A similar user option to input K f for the connecting pipe between 
the turbine is also included in the code. 
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6. WALL CONDUCTION MODELS 
6. 1 Intr od uet ion 

As described in sections 4.4 and 5.7, the heat transfer rates from the 
ga3 to the walls of the various system components depend on the instantaneous 
difference between the gas and the wall temperature. Estimates of these 
engine wall temperatures have been calculated in the past [42], [6] based on 
the following assumptions: 

(i) each heat transfer surface of interest is at a uniform surface 
temperature; i.e. surface temperature variations across a particular area 
are neglected. 

(ii) The heat transfer . urface temperatures are constant with time; i.e. 
cycle-periodic surface u ""oerature variations are not considered. 

(iii) Heat transfer by conduction through the walls is treated on a one- 
dimensional basis. 

The one-dimensional, uniform surface temperature model is certainly a 
simplification for surfaces with large temperature variations, such as piston 
bowls, or for surfaces with complicated heat transfer path lengths, such as 
the oyiinder liner. AI 30 , uniform surface temperatures are not adequate for 
detailed thermal stress calculations. For most surfaces, however, temperature 
varies much more rapidly ir. directions perpendicular to the surface, so that 
the above one-dimensional treatment is Justified. 

The assumption of constant wall temperatures over the engine operating 
cycle is reasonable for engines with high conductivity metal walls and forced 
convection Jacket cooling. Measurements by LeFeuvre [27] and Whitehouse [43] 
on conventional engines have suggested that cyclic surface temperature swings 
are fairly small, ranging from 5 to 15 Kelvins. However, for engine surfaces 
insulated with low-conductivity materials, such as ceramics, surface 
temperature swings are expected to be more critical [44], 

64 






7 


/ * 



Figure 6 3 how 3 a typical ceramic/metal composite wall structure. The 
heat transfer rate from t gas to the wall is a harmonic function of time, 
with a period of one engine cycle. Thi3 time-dependent boundary condition 
will set up periodic temperature waves that will propagate into the wall 
structure. Because of the relatively low thermal diffusivities of ceramics, 
these disturbances will only penetrate a small distance from the surface of 
the material, beyond which the temperature distribution is steady-state. 

These cyclic transients, should not be confused with engine start-up wall 
transients which die-away once steady-state engine operation is established. 
Cyclic transients are superposed on the steady-state conduction temperature 
profiles. Thus, their effect should be taken into account for the accurate 
prediction of maximum wall surface temperatures (important for lubrication 
considerations) and temperature distribution in the wall structures (important 
for material stress considerations). 

Assuming uniform material properties which do not vary with temperature, 
the total temperature (steady plus time-periodic) T(x,t), at any point x 
within the wall, and at any time t, will satisfy the heat conduction equation, 
i.e. 


where a is the thermal diffusivity of the material. The solution 
is developed in the following sections by decomposing the problem 
steady-state and time-periodic components. 


( 6 - 1 ) 

to eq. (6-1) 
into its 
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6.2 Steady- S tate Problem 


The steady-state heat transfer rate per unit inside surface area, Q , 

3 

through a composite structure is given by 

% * u(f w * V (6 " 2) 

where U * over-all coefficient of heat transfer based on inside surface area 
T * steady-state inside wall surface temperature 
T c * outside wall surface temperature, coolant temperature, or ambient 
temperature 

The cylinder head and the piston crown of the Insulated engine are 
modelled as flat composite walls, while the cylinder liner, manifolds, and 
connecting ducting are modelled as cylindrical composite walls. The 
appropriate expressions for the over-all heat transfer coefficient and steady- 
state temperature distribution for these two types of wall structure are 
summarized in sections 6.2.1 and 6.2.2. The inside wall steady-state surface 
temperature is deterrned through a heat balance between the cycle-averaged 
heat transfer rate from the gas to the walls and the steady-state wall heat 
conduction. This involves an iterative procedure which is described in 
section 6.2.3. 

6.2.1 Flat composite wall 

The overall coefficient of heat transfer, Up, for a flat composite wall 
is defined as 


1 

O' 

p 
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i-1 


L 
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i 

i 


where n *> numbar of flat layers 


- thickneoj of l layer 


(6-3) 
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k. = thermal conductivity of layer 

If the boundary condition on the outside wall surface is specified in 
terms of an outside heat transfer coefficient and a coolant or ambient 
temperature (versus a specified outside wall temperature), a modified overall 
heat transfer coefficient is defined as 


1 n L. 1 

— = y -i f — 

U k. h 

p 1=1 1 c 


(6-4) 


where h is the heat transfer coefficient to the outside, 
c 

The steady-state temperature distribution for the i^ 1 layer is obtained 
from the solution of Eq. (6-1) in its steady form, i.e. 


9 2 T 

■ 0 


(6-5) 


The boundary conditions are 


VWi 1 ' at x “ x i 


(6-6a) 


T s 3 V Q s’W’ 3t X = Vi 


(8-6b) 


th 

leading to the following temperature distribution within the i layer 

y*> ■ ty-i.,) * W (6 - 7) 

6.2.2 Cylindrical coiqposite wall 

The overall coefficient of neat transfer, U c , based on the inside surface 
area of a cylindrical wall structure is defined as 


-L , ? l n l r Ati 'S aI 

r U ^ k 
To i»1 i 


( 6 - 8 ) 
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approximate estimates of these temperatures are assumed. Based on these 
estimates, the instantaneous heat transfer rates, convective and radiative, to 
the combustion chamber walls are calculated throughout the cy^le. At the end 
of the cycle, a heat balance is performed between the cycle-averaged gas /wall 

heat transfer rate and the heat conducted through the walls of each component 

to compute a new surface temperature. These new temperatures are used in the 
next cycle iteration, until the calculation converges. Details of the heat 
balance follow. 

The instantaneous heat transfer rate to the combustion chamber walls is 
calculated from 

Q (t) ■ h(t)(T (t) - T ) (6-13) 

w gw 

where T w = inside wall temperature at cycle start 
T^(t) = instantaneous gas temperature 

h(t) = instantaneous heat transfer coefficient 

During combustion, h(t) is replaced by an effective linearized heat 
transfer coefficient h ft ^(t), given by Eq. (D-7), to take into account the 
effects of radiation on the total heat transfer rate. 

The cycle-averaged gas/wall heat transfer rate is given by 

/ h(t)(T (t) - T )dt 

»» — <<H ‘ I) 

where f denotes Integration over the complete cycle. Equation (6-1 A) can be 
rewritten as 
» 

Q w * h T g " V* (6-15) 

where the bar denotes cycle-averaged. 
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The steady-state heat transfer rate conducted through a composite wall 
structure was expressed by equation (6-2), involving an over all wall heat 
transfer coefficient. Combining equations (6-2) and (6-15), an updated wall 
temperature can be computed according to 

• h~T + UT 

T - - — 8 ( 6 - 16 ) 

h + U 

where the prime denotes the temperature to be used in the next cycle 
iteration. The calculation is repeated until the new wall temperature is 
within a certain percentage of its value at the end of the previous cycle. 


6.3 Time-Periodic Problem 
6.3.1 Formulation of finite difference scheme 

The time-periodic part, T p (x,t), of the temperature distribution within 
any parallel slab will satisfy the unsteady conduction equation, i.e., 


3 2 T 1 3T 

— E „ . __B 

3x 2 “ 3t 


(6-17) 


In order to solve this continuous partial differential equation, a finite 
difference technique will be used. The slab is modeled by a number of N 
discrete nodal points, x^. At each node, the finite-difference approximation 
to the governing equation provides an algebraic equation connecting the 
instantaneous nodal temperature to those at the surrounding nodes. 

The finite difference schemes which have been used to solve partial 
differential equations of the form of Eq. (6-17), can be grouped into two 
general categories: explicit and implicit [453. In explicit schemes, the 
instantaneous value of the variable at any node is given by the values of the 


♦ 
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variable at the node and its neighboring nodes calculated from the previous 
time step. These schemes are fairly simple, and computationally very 
efficient. However, they are stable only under certain conditions and are 
limited to uniformly-spaced discrete nodes. On the other hand, implicit 
schemes are unconditionally stable, and some can handle arbitrarily spaced 
discrete nodes. From the computational point of view, though, implicit 
schemes are very intensive. Since the instantaneous value of a nodal variable 
depends on the neighboring values of the variable at the current and previous 
time steps, these schemes, unfortunately, involve large matrix inversions. 

For the purposes of the current work, it is desired to calculate the 
temperature distribution within the walls of the engine combustion chamber in 
parallel with the engine simulation calculation. Thus, special effort was 
made to develop a numerical scheme that would be: 

(i) least demanding in computer time 

( i i) able to handle arbitrarily-spaced discrete nodes, so as to maximize 
the accuracy of the solution within the relatively thin penetration 
depth of the cyclic engine transients 
(iii) stable. 

To satisfy all the above requirements, an Euler explicit scheme was 
suitably modified, by mapping the arbitrarily-spaced discrete noues x, of the 
desired solution domain into corresponding uniformly-spaced nodes y^, 
according to the following transformation: 


a + bx ,, 

y - (6-i8) 

This transformation w 3 selected because of its following merits [H6]: 
(i) three of the nodes x^ can be mapped into three specified nodes y^. 
(!i) all other nodes x, are then mapped smoothly in between the three 
specified nodes y^. 
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(ill) the transformation is differentiable up to any order. 

Using (6-18), Eq. (6-17) nan take the following form for earh node x,: 


[ dx! ‘3y^3x JJ , a 9t , 


(6-19) 


3y approximating the derivative with respeet to y through finite differences, 
and recalling ( 6— 1 8 ) , Eq. (6-19) can be re-written as 

(~)(&) [( 2 } “ (S) (P> ] = 1 (P) (6-20) 

Ay 3 x l 3x 1+1/2 3y +1/2 3x i - 1 /2 9y i-1/2 “ 3t i 

where Ay is the uniform spacing between two adjacent nodes in the transformed 
coordinate system. 

Once again, by taking finite difference expression; for the derivatives 
with respect to y, and after some rearrangement, we get 

° DYDXM f o o o 

T i = T i ♦ [DYDXR(T i+1 - T { ) - DYDXL(T i - T^)] (6-21) 


where DYDXM = (|*) 


(6-22a) 


DYDXR - (|[) 


i+1/2 


(6-22b) 


DYDXL - (|[) 


1 - 1/2 


(6- 22c) 


Ay^ 

-f - : Courant number 
oAt 


(6-23) 


and the superscript 0 denotes temperature values at the previous time step. 
Note that the method is stable provided that [M7] 


CN 2 2 


( 6 - 2*0 


For time-step sizes of the order of one engine crank angle, and material 

_C .£ O 

diffusivities in the range of 10 5 to 10 m /s (most materials of interest), 
this condition will be satisfied for nodal spacings as small as 1/10 to 1/100 


of a mm. 
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6.3.2 App licat ion of f inite difference scheme 

The finite difference scheme was applied to calculate the periodic 
temperature distribution in the two parallel layer wall structure (piston or 
head), shown in Fig. 6. In general, under typical engine operating 
conditions, the cyclic transients do not penetrate into the wall structure for 
a distance greater than a few mm (and thus will not extend into the second 
layer). However, to maximize the flexibility of the code to analyze even thin 
coatings (of order of a mm), the code is set-up with a capability to calculate 
transients over the entire wall region. Obviously, more nodes 3 hould be 
placed in the first layer than in the second one. 

For the first layer’, the transformation (6-18) is introduced, so that 
x * - L 1 maps into y » 0 

x * 0 maps into y => 1 

x ■ - L 1 +5 maps into y * F 

with F ■ 0.4 to 0.5, meaning that 40 to 50 percent of the nodes are placed 

within a distance of the order of the penetration depth 6 in the first layer. 

These three conditions determine a, b, c in Eq. (6-18). Then, the temperature 

T . at any node i of the first layer, can be calculated from Eq. (6-21). 

' 9 i 

For the second layer, a relatively small number of uniformly-spaced nodes 
is sufficient. Then, with y - x, Eq. (6-21) for the temperature 7, , at any 

^9 J 

node j of the second layer, simplifies to the following standard form [47]: 



* (0M) \i * lj.1 

CN 


(6-25) 


Note again that for stability, CN i 2. (The limit CN » 2 would imply that the 
local node temperature has no effect on its future value). 

At the interface between the two layers (x-0), continuity of temperature 
and heat flux leads to the following boundary conditions: 
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(6-26) 






* 


T = T 
1.N 2.1 


whe-'e N is number of nodes In first layer, and 


, ,i)T> . ,3T. 

“ K 2 ( 3x.J n 

x~0- x*0* 


or, introduring transformation (6-13) in the first layer, and aproximating the 
derivatives through finite differences, we get 


3y 3T, M - MT 


♦ T. 


a.. 


' 3x 


Ay 


- k„ 


-3T ♦ HT - T 

AX- 


(6-27) 


where Ay was defined in Eq. (6-20) and Ax 0 is the uniform spacing between two 
adjacent nodes in the second layer. 

To validate the numerical 3oheme .and to optimize the number and spacing 
of nodal points, the method was applied to calculate the temperature 
distribution in a two layer composite, slab subject to a harmonically varying 
gas temperature at x * - ^ , and a constant ambient temperature at x » 

The exact solution was also obtained analytically for the same problem [10]. 
The approximate results were shown to be in excellent agreement with the exact 
-esults for a range of different materials ana wall structures [10]. 


6.4 Combined Solution 

Tne total temperature distribution in each layer can be obtained by 
superposing the steady-state solution given by Eq. (6-7), and the time- 
periodic solution, given by Eqns. (6-21) and (6-25). The total temperature 
T'x,t) must satisfy the full boundary conditions at the inside and outside 
wall surfaces. 

At the hot gas side (< « - l^), 


- k,(~) ♦ J - h (T - T ) 

1 * x _ _ h s g g w 


(6-28) 
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where Q is the steady-state heat flux given by (6-2). h is an effective 
s g 

linearized coefficient (convective and radiative) defined by (D-7), and 


T ♦ T. , 
w 1,1 


(6-29) 


Introducing finite-difference approximations, (6-28) can be rewritten as 

( 6 - 30 ) 


k (&) 

r3x J 


3T - 4T ♦ T 

3 1,1 1.2 


-L 


LIS 
2Ay 


\ l ! t q »= h (T - T ) 
s g g w 


Similar equations apply for the boundary condition at the ambient or 
coolant side (x - L 0 ), where either the total heat f jc or the total wall 
temperature can be prescribed. 

The calculation is staged as follows. First, the steady-state solution 
for the temperature distribution within each component wall and the average 
heat flux conducted through each wall are obtained by following the iteration 
procedure described in section 6.2.3, i.e. assuming no cyclic transients. 

Once the steady-state heat flux is obtained, the transient calculation is 
performed by applying the finite difference scheme and the total boundary 
conditions (6-28) and (6-30). 
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7. METHOD OF SOLUTION AND PROGRAM INPUTS AND OUTPUTS 
7 . 1 8a3tc Method of Solution 

The conservation equations of mass and energy for the contents of an open 
thermodynamic system are applied in turn to the master reciprocate.* cylinder, 
the intake manifold, and the series of exhaust manifold sections. Further, 
the individual submodels of the various system components and their 
thermodynamic and heat transfer processes are brought together to form a 
complete system model. The result is a set of simultaneous first-order 
ordinary differential equations. To perform predictive calculations with the 
cycle simulation, these equations mu3t be integrated simultaneously over the 
full operating cycle. Note, however, that some of the governing equations 
like the mass flow rate through the Intake or the exhaust valve, the non- 
dimens Lonal fuel burning rate, etc. apply only during parts of the cycle. 

Integration of the governing equations is performed numerically using a 
standardized code developed by Sharp inc and Cordon [48]. The method is based 
on a predictor-corrector technique that uses a modified divided difference 
form of the Adams Pace formulas. The code adjusts it3 order and step size 
internally to maximize efficiency and control the local error per unit step in 
a generalized sence. Detailed documentation of the integration routine is 
provided m the listing of the code. 

A fLow chart showing the overall structure of the entire system 
simulation is shown in Figure 7. After the initialization of the state 
variables, the program proceeds with the simultaneous integration of the 
governing system equations. The latter can be grouped into two major sub- 
sets: equations describing the thermodynamic processes in the master cylinder 

of the reciprocato 1 **; and equations associated with other components in the 
system l manifolds, turbocharger, etc.) that have Inherent multi-cylinder 
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Fig. 7: Flow chart of simulation program. 




capability. The main program determines which equations in each sub-set will 
oe passed to the integrator, coupler the two sub- sets and transfers 
information across them, and prints out results as the cycle proceeds. 

On the reciprocator model side, the main program passes to the Integrator 
different subroutines corresponding to the intake, compression, combustion and 
exhaust processes of the master cylinder, as the engine cycle proceeds. On 
the other components side, the main program supplies the Integrator with the 
same set of governing equations for the manifolds, turbocharger, etc. 
throughout the cycle. Appropriate utility routines, are called by the major 
routines to help in the evaluation of the necessary derivatives at eacn step. 
These include routines to calculate the steidy-state and cyclic periodic 
temperature profiles withii the various component wall structures, routines to 
calculate mass flow rates through valves and interpolate valve area tables 
(called by the intake and exhaust routines), routines to interpolate 
turboraachinery maps and predict the heat transfer and pressure losses in the 
ducting (called by the other components routine), and the<"nrdynamic and 
transport property routines (called by all major routJ>»..i). 

The mass flow rate and enthalpy flux profiles of the master cylinder, 
generated by the simulation during the intake anc. the exhaust processes, are 
stored in the main program. Using this information, the main program calls a 
subroutine that: 

i) generates the profiles of all the other cylinders as echoes of the 
master cylinder profiles; 

11) sums the Intake and exhaust profiles contributed by oach cylinder at 
every instant to give total reciprocator mass and enthalpy flows, etc. The 
resultant profiles are communicated from the main program to the other 
c exponents routine throughout the engine cycle. 


Approximate estimates of all state variables are assumed, initially. 

More than one iteration is generally required to model system operation under 
steady conditions. The integration continues until the system reaches a 
quasi-steady condition, defined as the condition in which the value of each 
state variable at a particular crank-angle is within a specified interval of 
its value at that crank-angle in the previous engine cycle. 

7.2 Program Inputs and Outputs 

The input parameters which must be specified for each cycle simulation 
calculation include the following groups: system operating mode, system 
operating conditions, system dimensions and design parameters, parameters for 
the wall conduction models, empirical parameters for the various simulation 
sub-models, initial conditions of the system and certain computational 
par- neters. 

From this information, the simulation program can predict the performance 
of the total engine system under a wide range of operating conditions. The 
output includes mean engine performance parameters, such as power, specific 
fuel consumption, mean effective pressure, thermal eff ‘ iency, etc. , as well 
as detailed information about the state of the total system as a function of 
crank-angle tin oughout each engine cycle. 

7.2.1 In puts 

I. System operating mode 

In order to maximize the flexibility of the code to handle different 
system configurations, and assess the effect of different modeling assumptions 
on system performance , the user i3 provided with a choice of several options: 

a. The system can be turbocharged and turbocompounded, or turbocharged 
onl>. 
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b. The steady-state inside wall surface temperatures of the piston, 
cylinder head, cylinder line*', manifolds, and turbine connecting pipe 
can be specified as an input, or calculated from a specification of 
the component wall structure. 

c. Additionally, the time-dependent temperature distribution in the 
pistcn and cylinder head can be computed using a one-dimensional 
unsteady finite difference model for the component wall. 

d. The ignition delay period can be predicted based on our ignition 
delay model, or can be specifed based on experimental data. 

e. Either the commonly used Annand's radiation model, or the flame 
radiation model developed here can be used to predict instantaneous 
radiative heat transfer rates. 

II. System operating conditions 

a. Reciprocator speed (RPM) 

b. Mass of fuel injected per cycle 

c. Injection timing 

d. Fuel parameters 

e. Steady-state inside wall surface temperatures of piston, cylinder 
head, cylinder liner, lanifold walls, and turbine connecting pipe 
(for specified wall temperature option only) 

f. Fower turbine gear ratio 

III. System dimensions and design parameters 

i. Reciprocator Parameters 

a. Number of cylinders 

b. Cylinder bore 


Cylinder stroke 


d. Connecting rod length 

e. Clearance volume 

f. Valve timings (crank angles at which the inta. e and 
exhaust valves open and close) 

g. Tabulated values for the effective valve open areas 
(including disciiarge coefficient effects) 

ii. Other Component Dimensions 

a. Intake and exhaust manifold dimensions 

b. Turbine connecting pipe dimensions 

iii. Intercooler Characteristics 

a. Coolant flow heat capacity 

b. Heat exchange surface area 

c. Over-all heat transfer coefficient 

iv. Turbomachinery Parameters 

a. Compressor, turbocharger turbine and power turbine maps 

b. Turbocharger rotational inertia 

c. Turbocharger rotational damping 

d. Powr turbine transmission efficiency 
IV. Wall conduction model parameters 

These parameters need to be specified when it is desired to predict the 
steady-state and transient (for pi3ton and cylinder head only) temperature 
distribution within the various system component walls. For each material 
layer of the piston, cylinder head, liner, manifold sections and turbine 
connecting pipe, the wall structure specifications include: 

a. Thickness 

b. Inside wall radius (cylindrical components only) 
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c. Thermal conductivity 

d. Thermal diff'-sivity (piston and head only) 

For each component, the boundary conditions on the outside wall surface 
must be specified, i.e.: 

a. Ambient (or coolant) temperature, and heat transfer 
coefficient from the outside wall surface to the ambient 
(or coolant), or 

b. specified W3ll temperature on outside surface 
finally, for the finite difference scheme applied to the unsteady 

conduction through the piston and cylinder wall, the following are required 

a. Number of nodes placed within each layer of a given component 

b. Fraction of nodes of first layer placed within the penetration 
depth of cyclic engine transients. 

V. Other sub -model empir ical paramete rs 

a. Ignition delay correlation constants; e.g., (4-12) 

b. Burning rate distribution constants; e.g., (4-7), (4-8), 

(4-9) 

c. Nu-Re number correlation constants, appropriate for the 
reciprocator cylinder, exhaust port, manifolds, and turbine 
connecting pipe; Equations (4-17), and (5-17) up to (5-21) 

d. Turbulent dissipation constant; Eq. (4-34) 

e. Annand's radiation model calibrating constant 

f. Friction model constants; «-.g., (4-48). 

IV. Initial conditions 

a. Cylinder pressure, temperature, and composition 

b. Intake manifold pressure, temperature, and composition 

c . Exhaust manifold pre33ure, temperature, and composition 
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d. Turbocharger speed 
VII. Ambient cond itions 

a. Intake temperature 

b. Intake pressure 

c. Final exhaust pressure 
VIII. Computational parameters 

a. Convergence margins for each state variable 

b. Error tolerances for integration of the differential equations 

c. Other parameters used in the integration algorithm "ODERT" 

(For a detailed description of these parameters, see the engine 
simulation code.) 

7.2.2 Outputs 

Four types of outputs are generated by the cycle simulation: 

I. Input echo 

A listing of all the input parameters, including some quantities derived 
directly from the given inputs (e.g. engine displacement and compression 
ratio). 

II. Major crank-angle by crank angle results 

At specified crank-angle intervals, the values of the following state 
variables are returned: 

a. Cylinder pressure, temperature, and average equivalence ratio 

b. Intake manifold pressure, temperature and average equivalence 
ratio 

c. Exhaust manifold pressure, temperature and average equivalence 
ratio 

In addition, the following other quantities are reported at the same 
intervals: 


84 



- -- k^r_' 


* 


4 




* 



« 



a. Heat transfer model results (such as heat transfer rates from 
gas to various component wall surfaces, instantaneous convective 
heat transfer coefficient, temperature profiles within various 
component wall structures) 

b. Turbulent flow model results (such as mean flow velocity, 
turbulent intensity, macroscale of turbulence) 

c. Code which monitors the performance of the integration routine. 
Integrating through the different processes for the master cylinder, the 
following quantities are reported during the corresponding process: 

i. Intake 

a. Mass flow through intake valve 

b. Mass flow through exhaust valve 

c. Velocity through intake valve 

d. Velocity through exhaust valve 

ii. Combustion 

a. Non-dimensional fuel burning rate 

b. Fuel burnt as a function of total fuel injected 

c. Flame radiation model results (such as radiant heat transfer, 
radiant temperature, adiabatic flame temperature, emissivity) 

iii. Exhaust 

a. Mass flow rate through exhaust valve 

b. Velocity through exhaust valve 
III. Integrated results and cycle performance 

After completion of an engine cycle, a summary of results obtained by 
intf rating some of the governing equations over the cycle is given. 

Integrated results include the following: 

a. Reciprocator thermal efficiency (gross indicated and brake) 





i J 
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b. Overall system brake thermal efficiency 

c. Specific fuel consumption (reciprocator and overall) 

d. Volumetric efficiency (based on manifold and ambient conditions) 

e. Gross indicated, punping, friction, and mean brake effective 
pressures 

f. Total heat loss (as a fractoin of the fuel energy input) 

g. Mean exhaust temperature 

h. Mass of air inducted per cycle 

i. Ignition delay period 

j. Total heat and work transferred during each process 

k. Results of an overall energy balance 
IV. Sub-model results 

After the overall cycle results are listed, detailed results for the main 
sub-models of the cycle simulation are given at specified crank angle 
intervals. These quantities include the following: 

a. Total engine intake and exhaust mass flow rates 

b. Compressor, turbocharger turbine and power turbine speed, mass 
flow rate, pressure ratio and efficiency 

c. Power turbine work transfer 

d. Pressures and temperatures at var.ous system locations 

e. Intake manifold, exhaust manifold, and turbine connecting pipe 
heat transfer data 

f. Intercooler performance data 
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APPENDIX A 

Thermodynamic Properties 

Our thermodynamic model assumes that the various open systems contain 
mixtures of air and combustion products throughout the total engine system. 

By utilizing the concept of the instantaneous average equivalence ratio 
defined in Section 4.1, the contents of <*ny open system can be represented as 
one continuous medium. Furthermore, assuming ideal gas behavior and 
thermodynamic equilibrium, the instantaneous gas properties can be determined 
from a knowledge of pressure, temperature and average equivalence ratio in the 
open system. 

When the terrperature of the cylinder contents is below 1000 K, they are 
treated as a homogeneous mixture of non-reacting lJeal gases, their properties 
being calculated using the procedure outlined below [49]: 

The hydrocarbon-air combustion reaction is written as: 

CH ♦ 0 2 ♦ * N 2 ♦ x 1 CO + x 2 H 2 0 ♦ x 3 C0 ♦ x^ ♦ x^ ♦ «N (A-1 ) 

v#iere ip » the molar N:0 ratio of the products, 
y - the molar 'i:C ratio of the fuel, 

$ * the average equivalence ratio, 

x. » moles of species i per mole of 0 ? reactant 

and e - 4/(4*y) (A-2) 

The quantities x, are determined by using the following assunptions: 

a) for lean mixtures ($ i 1) H ? can be neglected. 

b) for rich mixtures ( > 1 ) 0^ can be neglected. 

c) for rich mixtures, the gas water reaction 

C0 2 ♦ H 2 * GO ♦ H 2 0 

is in equilibrium with equilibrium constant K(T). 


2 Jk * hi- 



i. 

\\ 



V 


1 

1 


1 

1 P 

1 

A 






The solution for the x^ is shown in Table A-1 , where C Is obtained by 
solving equation (A-3) for its positive root. 

(1 - K)C 2 ♦ 2[ 1 - €♦ ♦ K(t - 1 ♦ e*)]C - 2M(* - 1) - 0 (A-3) 

The value of K(T) is obtained by curve fitting JANAF table data over the 
temperature range 400 to 3200 K and is given by 

ln(K(T)) - 2.743 - 1 .761 /t - 1.61 1/t 2 ♦ .2803/t 3 (A-4) 

vhere t - T/ICOO, and T is the temperature in Kelvins. 

If the grams of products per mole of 0^ reactant is expressed as 
M - (8e ♦ 4)<j> ♦ 32 ♦ 28* (A-5) 

the specific enthalpy h and the specific heats at constant pressure and 


constant composition, c and c respectively, can be expressed by the 

P * 


following relationships: 

1 6 U t J a i5 

h - „ E x. E (a. . - - -■* ♦ a.,) 

M i.i 1 j.i J 1 16 

c - n E x E (a t J ♦ ~*> 

p M l -1 *>1 lJ t 2 

1 6 3x t 4 j a 

” H ^ 3*' I " 't 5 * a i6 } " 


(A-6) 


(A-7) 


- E x E (a * a ) 

« 2 w 1 jit « j 1 16 


(A-8) 


The coefficients a,j are obtained by curve fitting JANAF table data to the 
above functional form. The values of a^ are given in Table A-2. The 

resultant c p Is In oal/g-K, while h and c^ are in kcal/g. 

Since the cylinder contents are being treated as a mixture of non- 
eactlng ideal gases, the density of the mixture is given by 
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TABLE A-1 

Burned Gas Compoc't^n under 1000 K 


i Spec ies 


x, (moles/mole 0 2 reactant) 


‘l 

4 > < 1 


<p > 1 


1 

2 

3 

4 

5 

6 


CO. 


h 2 o 

CO 


N 2 

Sun 


c<t> 

2(1-cH 

0 

0 

* 


e* - C 
2(1-e$) ♦ C 
C 

2(*-1) - C 
0 


(1-e)t ♦ 1 ♦ i|i (2-e)$ ♦ i|» 


Source: Ref. [49] 
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TABLE A-2 


oefficients for Polynomial Fit to Thermodynamic Properties 


Coefficients for 100 K < T < 500 K: 


i 

Species 

a il 

a i2 

a 13 

a i4 

a i5 

A * 

a i6 

2 

°o 2 

4.7373 

16.653 

-11.232 

2.8280 

.006767 

-93.75793 

3 

-2° 

7.8097 

-.20235 

3.4187 

-1.1790 

.001436 

-57.08004 

3 

co 

6.9739 

-.82383 

2.9420 

-1.1762 

.0004132 

-27.19597 

4 

»2 

6.991$ 

.16170 

-.21821 

.29682 

-.016252 

-.118189 

5 

°2 

6.2957 

2.3884 

-.031479 

-.32674 

.004359 

.103637 

6 

N 2 

7.0922 

-1.2958 

3.2069 

-1 .2022 

-.0003458 

-.013967 

Coefficients for 500 K 

< T < 6000 

K: 




i 

Species 

a i1 

a i2 

a i3 

a i4 

a i5 

a i6* 

1 

°°2 

11.940 

2.0886 

-.47029 

.037363 

-.58945 

-97.1418 

2 

h 2 o 

6.1391 

4.6078 

-.93560 

.066695 

.033580 

-56.62588 

3 

00 

7.0996 

1 .2760 

-.28775 

.022356 

-.15987 

-27.73464 

4 

H 2 

5.5557 

1 .7872 

-.28813 

.019515 

.16118 

.76498 

5 

°2 

7.8658 

.68837 

-.031944 

-.0026870 

-.20139 

-.893455 

6 

N 2 

6.8078 

1.4534 

-.32899 

.025610 

-.11895 

-.331835 


* picked to give enthalpy datum at 0 K 
Source: Ref. [49] 
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Air 



where 

R q - the universal gas constant (1.9869 oal/mole-K) 

and M, the average molecular weight of the mixture, is given by 
M - M/C (1 - c)<J> + 1 * i|») -j> < 1 

(A-10) 

M * M/( (2 - e)<1> ♦ i|») $ > 1 

Then, the partial derivatives of the density with respect to 
temperature, pressure, and equivalence ratio are given by 


lQ . . 9 
BT T 


(A-1 1 ) 


„ fi 

9p p 


(A-1 2) 


Bg Bg 3M g 3M ,, 15 , 

R'fiS 

When the temperature of the cylinder contents is above 1000 K, their 
properties are calculated with allowance for chemical dissociation, according 
to the calculation method described in [50]. This is an approximate method 
based on curve fitting data obtained from detailed thermochemical calculations 
[30] to a functional form obtained from a consideration of carbon air 
combustion. Although species concentrations within the burned gases are not 
calculated, the bulk thermodynamic properties needed for cycle analysis are 
accurately determined. 
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APPENDIX B 

Thermodynamic Data for Fuel Vapor 


Ir the computer engine simulation code, the fuel vapor enthalpy is 
modelled by a polynomial of the form 


H(T) 


t 2 t 3 

V + A 2 2 + A 3 3" 


♦ A 


4 4 



(B-1) 


Where t * T/1000, T is the temperature in Kelvins, and the units of H(T) are 
kcal/gmole. 

The value of Ag for a particular fuel depends on the datum at which the 
elements C(graphite), 0 2 and N 0 are assigned zero enthalpy. Using a 

datum of 298 K, the coefficients A were obtained for fuel C, n 0 „H, 0 , Q by 
curve fitting table data from Rossini et al [52] to the above functional 
form. The values of A^ for various hydrocarbon fuels, including #2 Diesel 
fuel, are given in Table B-1 . 

Our thermodynamic property computer codes, however, use a OK datum. To 
convert to a OK datum, a correction term Ag must be added to the enthalpy 
given by (E-l ) to account for the enthalpy difference of C and H 2 between 0 
and 298 K. Using data from the JANAF Tables [53], 


A 8 " n(H 0 - H °298 >c * I (H 0 - H V 


- 10.84 x 0.252 + x 2.024 


21 .636 kcal/gmole 
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APPENDIX C 


transport Properties 

The heat transfer correlations relate the heat transfer coefficient to 
the Reynolds and Prandtl numbers and the thermal conductivity. The 

calculation of the heat transfer rates will therefore require values for the / 

viscosity and the Prandtl number (from which the thermal conductivity can be 

Obtained). We have used the approximate correlations for the viscosity ana 

the Prandtl number of hydrocarbon-air combustion products developed by 

Mansouri and Hey wood [51]. 

The NASA equilibrium program [ 30 ] was used to compute the viscosity of 
hydrocarbon-air combustion products as a function of temperature, T, 
equivalence ratio, (|», and pressure p. It was shown that the viscosity of the 
combustion products was satisfactorily correlated by a power-law based on air 
viscosity data, corrected for the effect of equivalence ratio, i.e., 

Mprod Ckg/ms] - 3.3x10 -ir T 0,7 /(1 + 0.027*) 

(C-1) 

for 500 K i T i H000 K and 0 S <|> S H 
Note that the viscosity of the combustion products is independent of the 
pressure. 

The equilibrium Prandtl number of hydrocarbon-air combustion products was 
also calculated over the above ranges of temperature, pressure, and 
equivalence ratio. Using a second order polynomial of Y to curve fit the 
above data, it was shown that the following correlation for lean (<|> < 1) 

$ 

mixtures predicted values in good agreement (within 5 %) with the data, i.e., 

Pr - 0.05 ♦ H.2(Y - 1) - 6.7(Y - I) 2 (C-2a) 

for 500 K S T S A000 K and t S 1 
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For rich mixtures (4> > 1), a reasonable fit (les3 than 1 0% error) to the 
equilibrium Prandtl nunfoer values calculated with the NASA program was found 
to be the following: 

Pr - [0.05 + *I.2(Y - 1) - 6.7(Y * 1) 2 ]/[1 ♦ 0.015 x10‘ 6 (1>T) 2 ] 

for 2000 KST^ 3500 K and 1 < <p S 4 (C-2b) 


( 

T 
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APPENDIX D 


Linearization of Total Heat Transfer Rate at the Gas/Wall Interface 


The total instantaneous heat transfer rate, convective and radiative, 
per unit area, to the combustion chamber walls is given by 


Q - Q ♦ Q 
w c r 


where <2 - h(T - T ) 
o gw 


and ^ * k/lj* - tJ*) 

where the symbols used have the same meaning as in Section 6. A. 
The radiative heat transfer can be expressed as 

4 J| 4 4 

0 - k (T - T ) ♦ k (T - T ) 
r r r g r g w 


(D-1) 

(D-2) 

(D— 3) 


(D-*4) 


or alternatively as 


li ii 

(T* " T ) 


Cl » k -C- r --£-- (T . T ) + k (X 3 + x V + T T 2 + T 3 ) (T - T ) 
r r T - T g w r g gw gw w g w 

g w 

(D-5) 

Combining equations (D-1), (D-2) and (D-5), the total instantaneous heat 
transfer rate per unit area can be expressed in linear form as 


A.- w w 


(D-6) 


where h @ff is an effective linearized heat transfer coefficient defined as 


ii U 

T - T ) 

K rf - h + k (T 3 + T ^T ♦ T T 2 + T 3 * -“«*--) (D-7) 

err r g gw gw w t-t 

g w 
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